11institutetext: The University of Tokyo, ICRR, Kamioka Observatory and Kavli
IPMU, Kamioka, Gifu, Japan 22institutetext: Michigan State University,
Department of Physics and Astronomy, East Lansing, Michigan, U.S.A.

# The Neutrino Interaction Simulation Program Library

Yoshinari Hayato 11<EMAIL_ADDRESS>Luke Pickering 22

###### Abstract

is a neutrino-nucleus interaction simulation program library. It can be used
to simulate interactions for neutrinos with between 100 MeV and a few TeV of
energy. is also capable of simulating hadron interactions within a nucleus and
is used to model nucleon decay and hadron–nucleus interactions for particle
propagation in detector simulations. This article describes the range of
interactions modelled and how each is implemented.

##  1 Introduction

is primarily a neutrino–nucleus scattering simulation program library and
provides a complete model capable of predicting the observations for a wide
range of neutrino scattering experiments. is capable of simulating
neutrino–nucleon and coherent neutrino–nucleus interactions in a number of
reaction channels over a neutrino energy range from 100 MeV to a few TeV.
Additionally, incorporates initial- and final state nuclear effects for
interactions with nuclei from boron to lead. One of the most important nuclear
effects is the re-scattering of hadrons, which are produced in the primary
neutrino–nucleon interaction, as they propagate out of the nuclear medium.
This re-scattering can result in hadron absorption, extra hadron production or
knock-out, or distortion of the nuclear-leaving particle kinematic spectra.
The hadron re-scattering model has also been used to simulate low-energy
pion–nucleus scattering both to tune the model to extant data Elder:2019  and
to simulate pion propagation in neutrino-scattering experimental simulations.
Finally, can also simulate various nucleon decay channels to support
experimental searches for the process.

has a long rich history, originally developed in the 1980s as a tool to study
atmospheric neutrinos and nucleon decay in the Kamiokande experiment
Kamiokande:1986iox , and some of the original FORTRAN77 code is still in use.
continues to be predominantly developed and maintained by members of the
Kamiokande series of experiments (Super-Kamiokande, T2K, Hyper-Kamiokande) and
many source files contain comments messages from the numerous physicists who
have contributed to the simulation over the past 35 years—including those
working on the Nobel prize-winning Super-Kamiokande (SK) analysis Super-
Kamiokande:1998kpq . Recent development has targeted the improvements most-
needed for precise and robust analyses of SK and T2K neutrino oscillation data
and neutrino cross-section measurements. At T2K energies, Charged-current
quasi-elastic interactions dominate. It has become clear over the last decade
that such interactions can only be precisely predicted by incorporating
detailed models of relevant nuclear dynamics. For the multi-GeV samples used
in SK analyses, shallow and deep inelastic scattering channels are critical
for modelling the expected rate of multi-ring events seen in the detector. The
transition region between resonance excitation and deep inelastic scattering
has proven particularly difficult to model well. Because of the in-house
nature of development and analysis usage, it is not yet open source. We do not
yet have the resources to migrate to an open source model, but, access to the
code and usage instructions are available upon request.

##  2 Motivating the use of an Interaction Simulation

The primary goal of atmospheric and long-baseline neutrino-scattering
experiments is to study neutrino oscillation, which occurs as a function of
neutrino energy and flavor. As neutrinos are neutral particles, their
properties (including their energy) can only be inferred from observable
secondary particles produced when they interact with matter in our particle
detectors. Interaction simulations are used to predict the probability of
neutrinos of a given flavor and energy to interact with a given target, the
observable secondary particle spectra produced in the neutrino–nucleus
interaction, and event selection efficiencies and purities. The common use of
nuclear targets in neutrino scattering experiments further complicates the
problem as secondary hadrons produced in the _primary_ neutrino interaction
can be absorbed or lose energy before leaving the target nucleus, obfuscating
the details of the primary interaction. As a result, interaction simulations
are a critical tool in the analysis of neutrino oscillation data. The neutrino
energy range modeled most carefully by is between 0.1 and 10 GeV. This range
is motivated by the typical energies for oscillation features expected in SK
and T2K data, as seen in Fig. 1 and Fig. 2.

Figure 1: Oscillation probabilities for atmospheric $\nu_{\mu}$ to $\nu_{e}$
(top) and $\bar{\nu}_{\mu}$ to $\bar{\nu}_{e}$ (bottom) as a function of
neutrino energy and the cosine of the zenith angle (acting as a proxy for
propagation distance through the Earth) assuming a normal neutrino mass
hierarchy. Matter effects in the Earth produce distortions in oscillations of
$\nu_{\mu}$ between two and ten GeV, which are not seen for $\bar{\nu}_{\mu}$
oscillations. Figures reproduced from Ref. Super-Kamiokande:2017yvm  The
oscillation parameters used are: $\Delta m_{32}^{2}=2.5\times
10^{-3}\mbox{eV}^{2}$, $\mbox{sin}^{2}\theta_{23}=0.5$,
$\mbox{sin}^{2}\theta_{13}=0.0219$, and $\delta_{CP}=0$.

Figure 2: The T2K $\nu_{\mu}$ flux T2K:2015sqm ; T2K:2012nrr  (energy spectral
shape) and overlaid oscillation probabilities for the $\nu_{\mu}$ to
$\nu_{\mu}$ (top) and $\nu_{\mu}$ to $\nu_{e}$ and $\nu_{\mu}$ to $\nu_{e}$
(bottom) as a function of neutrino energy. The black (blue) curves correspond
to minimal (maximal) values of the CP-violating phase $\delta_{\textsc{cp}}$.
The oscillation parameters used are: $\Delta m^{2}_{32}=2.5\times
10^{-3}\rm{eV^{2}}$, $\sin^{2}(\theta_{23})=0.5$,
$\sin^{2}(2\theta_{13})=0.1$.

To test oscillation hypotheses, observable distributions that correlate
strongly with the neutrino energy distribution are sampled and compared to
simulated predictions to extract oscillation parameter constraints. An example
of such an observable is

$E^{\textrm{Rec}}_{\textsc{qe}}\left(p_{\ell},\theta_{\ell}\right)=\frac{2M_{N,i}E_{\ell}-M_{\ell}^{2}+M_{N,f}^{2}-M_{N,i}^{2}}{2\left(M_{N,i}-E_{\ell}+p_{\ell}\cos{\theta_{\ell}}\right)},$
(1)

where $M_{N,i}$, $M_{N,f}$, and $M_{\ell}$ are the mass of the initial-state
nucleon, final-state nucleon, and final-state charged lepton respectively;
$E_{\ell}$, $p_{\ell}$, and $\theta_{\ell}$ are the energy, three-momentum,
and angle of the final-state charged-lepton respectively.
$E^{\textrm{Rec}}_{\textsc{qe}}$ is an unbiased neutrino energy estimator for
the charged-current quasi-elastic (CCQE) reaction off a free neutron,
$\nu_{\ell}+n\rightarrow\ell^{-}+p$, or a free proton,
$\bar{\nu_{\ell}}+p\rightarrow\ell^{+}+n$, which is only a function of the
kinematics of the final-state charged lepton. All operating neutrino-
scattering experiments use nuclear targets, where CCQE reactions occur
predominantly with bound nucleons. For such interactions,
$E^{\textrm{Rec}}_{\textsc{qe}}$ is biased by the potential energy associated
with the nucleon binding and is smeared by the Fermi motion of the nucleon
confinement111This effect can be seen in the offset between the initial
neutrino energy arrow and the peak of the $E^{\textrm{Rec}}_{\textsc{qe}}$
distribution in Fig. 3 (middle).. For interactions that exhibit strong
nuclear-effects or produce extra hadrons, here referred to as _CC non-QE_ ,
$E^{\textrm{Rec}}_{\textsc{qe}}$ is in general a significant underestimate of
the true neutrino energy. Fig. 3 (top) shows the prediction of the oscillated
rate of charged-current neutrino–oxygen interactions that produce no
observable pions, a CC0$\pi$ event _topology_. The CCQE and CC non-QE
components are shown separately to highlight how the shape of their
contributions differ between the $E_{\nu}$ and
$E^{\textrm{Rec}}_{\textsc{qe}}$ projections. For non-QE interactions, the
shape of the oscillation is largely smeared away. The CC non-QE component of
the CC0$\pi$ topology can be further broken down into interactions that either
do not produce secondary pions, or do produce secondary pions that are
subsequently re-absorbed before leaving the struck nucleus. Fig. 3 (middle and
bottom) illustrates the -predicted evolution of
$E^{\textrm{Rec}}_{\textsc{qe}}$ for a relevant range of discrete neutrino
energies. It can be seen that $E^{\textrm{Rec}}_{\textsc{qe}}$ exhibits
significant reconstructed energy feed down above about 1 GeV, where secondary
meson production is common. For higher energies, calorimetric estimators that
account for hadronic energy are more accurate. However, current neutrino
interaction models are often only predictive in lepton kinematics and
significant uncertainties should be assigned to the predictions of hadronic
particle spectra, especially for baryons and heavy mesons.

To accurately interpret observable distributions, interaction simulations are
relied upon to predict the rate and observable projections for neutrino
interactions over a range of energies and for a number of different target
nuclei as well as the migration of these primary interactions into observable
topologies. The next section details the nuclear dynamics, primary neutrino
interaction, and hadronic re-scattering physics models implemented in .

Figure 3: (top) The true and reconstructed neutrino energy for charged-current
interactions producing no observable pion in a typical oscillated T2K
$\nu_{\mu}$muon neutrino flux at SK T2K:2015sqm ; T2K:2012nrr . The
reconstructed energy feed down for non-Quasi-elastic interactions smears away
the oscillation signature that can be clearly seen in the true energy
spectrum. (middle) The reconstructed energy feed down for four mono-energetic
muon neutrino fluxes. The colored arrows show the true energy of the
corresponding reconstructed energy distribution. (bottom) The 3 GeV
distribution is separated by number of secondary mesons produced.

The most important proton decay channels, predicted by various grand
unification theories, produce a lepton and a meson, _e.g._ $p\rightarrow
e^{+}\pi^{0}$, $p\rightarrow\mu^{+}\pi^{0}$, or $p\rightarrow\bar{\nu}K^{+}$.
Neutrinos also often produce a lepton and one or more mesons when interacting
with a nucleon. There are plenty of atmospheric neutrinos, which have
sufficient energy to undergo such interactions, providing a significant
background for proton decay searches. The majority of neutrino interactions
also produce a baryon in addition to a lepton and a meson, but these may not
be visible in detectors such as SK. The typical predicted momentum of mesons
from both proton decays and atmospheric neutrino interactions are below
several hundred MeV/c and thus, the probabilities for them to scatter before
leaving the nucleus is high. Accurately predicting neutrino-induced
backgrounds and hadronic re-scattering within the nucleus is critical for the
sensitivity of proton decay searches.

##  3 The Physics in

###  3.1 Simulating an Interaction

In general, factorizes the simulation of an interaction of a neutrino with
flavour, $\ell$, and energy, $E_{\nu}$, into four discrete steps. First, a
specific interaction channel is chosen randomly with probability,
$P=\sigma_{T}^{i}\left(E_{\nu_{\ell}}\right)/\sigma_{T}^{\textrm{tot}}\left(E_{\nu_{\ell}}\right)$,
where $\sigma_{T}^{\textrm{tot}}\left(E_{\nu_{\ell}}\right)$ is the total
cross section and $\sigma_{T}^{i}\left(E_{\nu_{\ell}}\right)$ is the cross
section for the specific target nuclei, $T$, and channel, $i$, where $i$ is an
integer that identifies the interaction process and is defined in Table 1
(charged current) and Table 2 (neutral current). For neutrino–nucleon
interaction channels, the nuclear-target cross section is usually constructed
as $\sigma_{T}^{i}=Z\sigma_{p}^{i}+(A-Z)\sigma_{n}^{i}$, where $A$ and $Z$ are
the nucleon number and the proton number of the target nuclei and
$\sigma_{p}^{i}$ and $\sigma_{n}^{i}$ are the bound proton and bound neutron
cross sections. For historical reasons, free protons can be added to nuclear
targets to build simple molecular targets such as $\textrm{H}_{2}$O and CH.
Fig. 4 shows the water-target cross-section predictions separated into classes
of interaction channel.

Table 1: Neutrino charged-current scattering reactions simulated by , where p refers to a proton, n to a neutron, N to a neutron or proton, ${}^{A}\mathbf{X}$ to an entire nucleus, and $W$ to the invariant mass of the final state hadronic system. The $\nu_{\mu}$ and $\bar{\nu_{\mu}}$ cross sections per nucleon for each process at 600 MeV and 10 GeV are included, where B.T. stands for Below Threshold and highlights kinematically disallowed channels for 600 MeV neutrinos. For reference, the reaction identifier enumeration is included. CCQE (1p1h) cross-sections are for bound nucleon in oxygen and calculated using the spectral function model. Coherent pion production cross-sections are also for oxygen. |  |  $\sigma^{\nu_{\mu}}/10^{-38}\textrm{cm}^{2}$ /N |
---|---|---|---
Channel Name | Reaction | 0.6 GeV | 10 GeV | Enum
CCQE (1p1h) | $\nu_{\ell}+\textrm{n}\rightarrow\ell^{-}+\textrm{p}$ | 0.76 | 0.95 |  -11
| $\bar{\nu}_{\ell}+\textrm{p}\rightarrow\ell^{+}+\textrm{n}$ | 0.20 | 0.85 | -11
2p2h | $\nu_{\ell}+\textrm{N}\textrm{n}\rightarrow\ell^{-}+\textrm{N}\textrm{p}$ | 0.03 | 0.08 |  -21
| $\bar{\nu}_{\ell}+\textrm{N}\textrm{p}\rightarrow\ell^{+}+\textrm{N}\textrm{n}$ | 0.01 | 0.08 | -21
CCRes$1\pi^{+}$ | $\nu_{\ell}+\textrm{p}\rightarrow\ell^{-}+\textrm{p}\pi^{+}$ | 0.15 | 0.77 |  -11
| $\nu_{\ell}+\textrm{n}\rightarrow\ell^{-}+\textrm{n}\pi^{+}$ | 0.03 | 0.52 |  -13
CCRes$1\pi^{0}$ | $\nu_{\ell}+\textrm{n}\rightarrow\ell^{-}+\textrm{p}\pi^{0}$ | 0.04 | 0.39 |  -12
| $\bar{\nu}_{\ell}+\textrm{p}\rightarrow\ell^{+}+\textrm{n}\pi^{0}$ | 0.01 | 0.31 | -12
CCRes$1\pi^{-}$ | $\bar{\nu}_{\ell}+\textrm{n}\rightarrow\ell^{+}+\textrm{n}\pi^{-}$ | 0.02 | 0.63 | -11
| $\bar{\nu}_{\ell}+\textrm{p}\rightarrow\ell^{+}+\textrm{p}\pi^{-}$ | $5\cdot{}10^{-3}$ | 0.41 | -13
CCDif$1\pi^{+}$ | $\nu_{\ell}+\textrm{p}\rightarrow\ell^{-}+\textrm{p}\pi^{+}$ | B.T. | 0.03 |  -15
CCDif$1\pi^{-}$ | $\bar{\nu}_{\ell}+\textrm{p}\rightarrow\ell^{+}+\textrm{p}\pi^{-}$ | B.T. | 0.03 | -15
CCCoh$1\pi^{+}$ | $\nu_{\ell}+{}^{A}\mathbf{X}\rightarrow\ell^{-}+{}^{A}\mathbf{X}+\pi^{+}$ | $1\cdot{}10^{-3}$ | 0.04 |  -16
CCCoh$1\pi^{-}$ | $\bar{\nu}_{\ell}+{}^{A}\mathbf{X}\rightarrow\ell^{+}+{}^{A}\mathbf{X}+\pi^{-}$ | $1\cdot{}10^{-3}$ | 0.04 | -16
CCRes1$\gamma$ | $\nu_{\ell}+\textrm{n}\rightarrow\ell^{-}+\textrm{p}\gamma$ | $2\cdot{}10^{-4}$ | $2\cdot{}10^{-3}$ |  -17
| $\bar{\nu}_{\ell}+\textrm{p}\rightarrow\ell^{+}+\textrm{n}\gamma$ | $3\cdot{}10^{-5}$ | $1\cdot{}10^{-3}$ | -17
CCN$\pi$ | $\nu_{\ell}+\textrm{N}\rightarrow\ell^{-}+\textrm{N}^{\prime}+x\pi$ | B.T. | 0.85 |  -21
| $\bar{\nu}_{\ell}+\textrm{N}\rightarrow\ell^{+}+\textrm{N}^{\prime}+x\pi$ | B.T. | 0.62 | -21
| where $x>1$ |  |  |
| and $1.3<W<2.0\leavevmode\nobreak\ \textrm{GeV}$ |  |  |
CCRes$1\eta^{0}$ | $\nu_{\ell}+\textrm{n}\rightarrow\ell^{-}+\textrm{p}\eta^{0}$ | B.T. | 0.19 |  -22
| $\bar{\nu}_{\ell}+\textrm{p}\rightarrow\ell^{+}+\textrm{n}\eta^{0}$ | B.T. | 0.14 | -22
CCRes$1K^{0}$ | $\nu_{\ell}+\textrm{n}\rightarrow\ell^{-}+\Lambda+K^{+}$ | B.T. | 0.06 |  -23
CCRes$1K^{+}$ | $\bar{\nu}_{\ell}+\textrm{p}\rightarrow\ell^{+}+\Lambda+K^{0}$ | B.T. | 0.03 | -23
CCDIS | $\nu_{\ell}+\textrm{N}\rightarrow\ell^{-}+\textrm{N}^{\prime}+x\pi$ | B.T. | 4.53 |  -26
| $\bar{\nu}_{\ell}+\textrm{N}\rightarrow\ell^{+}+\textrm{N}^{\prime}+x\pi$ | B.T. | 1.30 | -26
| where $W>2.0\leavevmode\nobreak\ \textrm{GeV}$ |  |  |
Table 2: Neutrino neutral-current scattering reactions simulated by , where p refers to a proton, n to a neutron, N to a neutron or proton, ${}^{A}\mathbf{X}$ to an entire nucleus, and $W$ to the invariant mass of the final state hadronic system. The $\nu_{\mu}$ and $\bar{\nu_{\mu}}$ cross sections per nucleon for each process at 600 MeV and 10 GeV are included, where B.T. stands for Below Threshold and highlights kinematically disallowed channels for 600 MeV neutrinos. For reference, the reaction identifier enumeration is included. NCEL (1p1h) cross-sections are for bound nucleon in oxygen and calculated using the spectral function model. Coherent pion production cross-sections are also for oxygen. |  |  $\sigma^{\nu_{\mu}}/10^{-38}\textrm{cm}^{2}$ /N |
---|---|---|---
Channel Name | Reaction | 0.6 GeV | 10 GeV | Enum
NCRes$1\pi^{0}$ | $\nu_{\ell}+\textrm{n}\rightarrow\nu_{\ell}+\textrm{n}\pi^{0}$ | 0.03 | 0.15 |  -31
|  | $9\cdot{}10^{-3}$ | 0.13 | -31
| $\nu_{\ell}+\textrm{p}\rightarrow\nu_{\ell}+\textrm{p}\pi^{0}$ | 0.03 | 0.15 |  -32
|  | $9\cdot{}10^{-3}$ | 0.12 | -32
NCRes$1\pi^{0}$ | $\nu_{\ell}+\textrm{n}\rightarrow\nu_{\ell}+\textrm{p}\pi^{-}$ | 0.02 | 0.13 |  -33
|  | $5\cdot{}10^{-3}$ | 0.11 | -33
NCRes$1\pi^{+}$ | $\nu_{\ell}+\textrm{p}\rightarrow\nu_{\ell}+\textrm{n}\pi^{+}$ | 0.02 | 0.12 |  -34
|  | $5\cdot{}10^{-3}$ | 0.10 | -34
NCDif$1\pi^{0}$ | $\nu_{\ell}+\textrm{p}\rightarrow\nu_{\ell}+\textrm{p}\pi^{0}$ | B.T. | 0.01 |  -35
|  | B.T. | 0.01 | -35
NCCoh$1\pi^{0}$ | $\nu_{\ell}+{}^{A}\mathbf{X}\rightarrow\nu_{\ell}+{}^{A}\mathbf{X}+\pi^{+}$ | $1\cdot{}10^{-3}$ | 0.02 |  -36
|  | $1\cdot{}10^{-3}$ | 0.02 | -36
NCRes1$\gamma$ | $\nu_{\ell}+\textrm{n}\rightarrow\nu_{\ell}+\textrm{n}\gamma$ | $2\cdot{}10^{-4}$ | $1\cdot{}10^{-3}$ |  -38
|  | $6\cdot{}10^{-5}$ | $8\cdot{}10^{-4}$ | -38
| $\nu_{\ell}+\textrm{p}\rightarrow\nu_{\ell}+\textrm{p}\gamma$ | $2\cdot{}10^{-4}$ | $9\cdot{}10^{-4}$ |  -39
|  | $6\cdot{}10^{-5}$ | $8\cdot{}10^{-4}$ | -39
NCN$\pi$ | $\nu_{\ell}+\textrm{N}\rightarrow\nu_{\ell}+\textrm{N}^{\prime}+x\pi$ | $2\cdot{}10^{-4}$ | 0.25 |  -41
|  | $3\cdot{}10^{-5}$ | 0.23 | -41
| where $x>1$ |  |  |
| and $1.3<W<2.0\leavevmode\nobreak\ \textrm{GeV}$ |  |  |
NCRes1$\eta^{0}$ | $\nu_{\ell}+\textrm{n}\rightarrow\nu_{\ell}+\textrm{n}\eta^{0}$ | B.T. | 0.04 |  -42
|  | B.T. | 0.03 | -42
| $\nu_{\ell}+\textrm{p}\rightarrow\nu_{\ell}+\textrm{p}\eta^{0}$ | B.T. | 0.03 |  -43
|  | B.T. | 0.02 | -43
NCRes1$K^{0}$ | $\nu_{\ell}+\textrm{n}\rightarrow\nu_{\ell}+\Lambda+K^{0}$ | B.T. | 0.01 |  -44
|  | B.T. | $7\cdot{}10^{-3}$ | -44
NCRes1$K^{+}$ | $\nu_{\ell}+\textrm{p}\rightarrow\nu_{\ell}+\Lambda+K^{+}$ | B.T. | $9\cdot{}10^{-3}$ |  -45
|  | B.T. | $7\cdot{}10^{-3}$ | -45
NCDIS | $\nu_{\ell}+\textrm{N}\rightarrow\nu_{\ell}+\textrm{N}^{\prime}+x\pi$ | B.T. | 1.37 |  -46
|  | B.T. | 0.49 | -45
| where $W>2.0\leavevmode\nobreak\ \textrm{GeV}$ |  |  |
NCEL (1p1h) | $\nu_{\ell}+\textrm{n}\rightarrow\nu_{\ell}+\textrm{n}$ | 0.13 | 0.16 |  -51
|  | 0.05 | 0.16 | -51
| $\nu_{\ell}+\textrm{p}\rightarrow\nu_{\ell}+\textrm{p}$ | 0.16 | 0.21 |  -52
|  | 0.06 | 0.19 | -52

Figure 4: The -predicted muon neutrino–water cross sections overlaid on the
T2K muon neutrino flux T2K:2012nrr , with a typical oscillation (top), and the
T2K and upward atmospheric muon neutrino fluxes Honda:2015fha  multiplied by
the charged-current inclusive total cross section (bottom). The flux
multiplied by the cross section is proportional to the expected interaction
rate. Above 4 GeV, the expected number of interactions in SK arising from the
T2K beam falls significantly faster than from atmospheric neutrinos. _n.b._
The cross sections presented in the top pane are divided by the neutrino
energy, whereas in the bottom pane, they are not. This is to emphasise the
saturation of the interaction channels associated with lower four-momentum
transfer at SK energies and the sharp turn-on seen over T2K flux distribution.

Second, the primary neutrino interaction, or _hard scatter_ , is simulated.
For the majority of channels, this step involves choosing a bound nucleon from
an initial-state nuclear model, then choosing interaction kinematics according
to the specific interaction model, and finally choosing any remaining particle
kinematics not specified by the model. This step is performed under the
impulse approximation PhysRevD.72.053005 , which treats the target bound
nucleon and the remnant nucleus as evolving independently during and after the
hard scatter. This further factorizes the simulation as, to first order, the
sampling of the nuclear model does not depend on the interaction kinematics
chosen.

For the coherent pion-production channels (Enum 16 and 36), the interaction
occurs coherently between the neutrino and the target nucleus and as a result
no bound nucleon target is chosen and this is considered the final step of the
simulation. For other channels, the final state hadrons are then passed on to
the third step, the nucleon and meson intra-nuclear re-scattering simulation,
where hadrons can elastically scatter, exchange charge with a nucleon in the
nucleus, or be produced or absorbed as they are stepped out of the nuclear
medium.

Finally, for oxygen targets only, the final state nuclear remnant can be left
in an excited state after the interaction and a number of nuclear de-
excitations, producing low energy photons
($\mathcal{O}\left(1-10\right)\,\textrm{MeV}$), are modeled following Ref.
PhysRevC.48.1442 . Careful treatment of the de-excitation oxygen is important
for precisely simulating interactions in the sensitive SK detector.

For the majority of particles produced in the hard scatter and subsequent re-
scattering, stores their properties in an event vector file that can be used
as input to further experiment simulation processes. The only exceptions are
tau and omega particles, which are decayed during the simulation by
TAUOLAJadach:1993hs  and a custom omega decay simulation code, respectively.
The final states of such decays are written to the event vector.

###  3.2 The Hard Scatter

This section details the neutrino scattering physics models implemented in ,
roughly ordered by the degree of inelasticity of the interaction.

####  3.2.1 Quasi-Elastic Scattering

The cross-section of neutrino–nucleon charged-current quasi-elastic (CCQE)
interaction was formalized by Llewellyn Smith QECRSFN . In some of the
literature, an equivalent channel is referred to as 1p1h, for one-particle
one-hole. Three different nuclear models are implemented in for simulating
CCQE interactions, the relativistic global Fermi-gas (GRFG) model, the local
Fermi-gas (LFG) model and the spectral function (SF) model.

The implementation of the GRFG coss-section follows the prescription by Smith
and Moniz QECRSFRMGAS . The LFG implementation uses the model by Nieves  _et
al._ Gran:2013kda  and includes an updated removal energy treatment
implemented by Bourguille  _et al._ Bourguille:2021:04 . Simple Fermi-gas
models tend to over-predict the cross section for forward going leptons, as a
result this model takes into account long and short range correlation of
nucleons using the random phase approximation, suppressing the cross-section
for low 4-momentum transfer Nieves:2004wx . The SF uses the spectral function
by Benhar  _et al._ Benhar:1994hw  and the implementation is based on the one
in NuWro NuWro2006  with additional improvements by Furmanski
Furmanski:2015knr .

These three nuclear models differ in their treatment of the bound nucleon
momentum and removal energy distributions and whether there are correlations
between them. Fig. 5 shows the projections of interactions simulated with each
of the three models into missing momentum ($p_{\textrm{miss}}$) and missing
energy ($E_{\textrm{miss}}$), which are observable quantities for analogous
measurements of electron–nuclear scattering. For the
$\nu_{\ell}+{}^{16}\textrm{O}\rightarrow\ell^{-}+p+{}^{15}\textrm{O}$
interaction, $p_{\textrm{miss}}$ and $E_{\textrm{miss}}$ are defined as

$\displaystyle p_{\textrm{miss}}$
$\displaystyle=\left|\vec{p}_{\nu}-\vec{p}_{\ell}-\vec{p}_{p}\right|,$ (2)
$\displaystyle T_{{}^{15}\textrm{O}}$
$\displaystyle=\sqrt{p_{\textrm{miss}}^{2}+M_{{}^{15}\textrm{O}}^{2}}-M_{{}^{15}\textrm{O}},$
(3) $\displaystyle E_{\textrm{miss}}$
$\displaystyle=E_{\nu}+M_{n}-E_{\ell}-E_{p}-T_{{}^{15}\textrm{O}},$ (4)

where $T_{{}^{15}\textrm{O}}$ is the reconstructed kinetic energy and
$M_{{}^{15}\textrm{O}}$ the ground-state mass of the nuclear remnant, which in
this case is an unstable isotope but is long-lived on the timescale of the
impulse approximation. These quantities are of interest because
$E_{\textrm{miss}}$ is approximately the energy lost to the nuclear response
during the interaction and $p_{\textrm{miss}}$ is the momentum of the struck
nucleon in the lab frame. It is clear from Fig. 5 that the different nuclear
models make different predictions about their distribution and correlation,
notably the ${}^{16}\textrm{O}$ shell structure is visible in the SF
predictions. The SF model is tuned to exclusive electron–nuclear scattering
data and is only available for a subset of target nuclei, and in is only
implemented for ${}^{12}\textrm{C}$, ${}^{16}\textrm{O}$, and
${}^{56}\textrm{Fe}$, the most important nuclei for SK and T2K analyses. For
interactions with other nuclei, such as ${}^{40}\textrm{Ar}$, only the GRFG or
the LFG nuclear models are implemented.

Figure 5: (Top) The reconstructed missing energy and missing momentum
distributions for the three initial-state nuclear models implemented for the
CCQE channel. (Bottom) Three nucleon axial form factor models and their
associated uncertainties derived in Ref. StowellThesis  (Dipole and
3-component) and Ref. PhysRevD.93.113015  (Z-expansion).

Beyond the choice of nuclear models, the vector and axial-vector nucleon form
factors control the strength and shape of quasi-elastic interactions. There
are two vector form factors of nucleon implemented in , the simple dipole form
factor and BBBA05 Bradford:2006yz . By default, BBBA05 is used as this form
factor was developed to reproduce experimental electron-scattering data.

There are four axial nucleon form factor models implemented, the dipole form
factor, BBBA07 Bodek:2007ym , the Z-expansion model Bhattacharya:2011  and the
3-component model StowellThesis . The 3-component fit model is inspired by the
2-component fit model Adamuscin:2007fk  and was created to provide additional
shape freedom by expanding the 2-component model, which quickly decays with
four-momentum transfer squared, $Q^{2}$. This model has the freedom to vary
the gradient of the form factor at low-$Q^{2}$ leaving the free parameters in
the 2-component model to set the shape at higher momentum transfer. The
3-component model is able to be continuously varied between the shape of both
the 2-component and the simple dipole model. Fig. 5 presents the shape of
three of the axial form factor models with associated uncertainties derived
fits to hydrogen and deuterium bubble chamber data. The dipole and 3-component
model uncertainties are reproduced from Stowell StowellThesis , and the
Z-expansion model was fit by Meyer  _et. al._ PhysRevD.93.113015 .

For neutral current elastic (NCEL) scattering, the treatment and available
model components are equivalent to CCQE, except that an LFG initial-state
model is not implemented for NCEL.

####  3.2.2 Charged Current Multi-Nucleon Scattering

A number of accelerator-based neutrino oscillation experiments, K2K, MINOS,
and MiniBooNE, started taking high-statistics neutrino–nuclear scattering data
in the 2000s. These experiments found that the number of the observed CCQE-
like (equivalent to the CC0$\pi$ topology introduced earlier) events were a
few tens of percent larger than predicted by the models, but with a relative
deficit of very forward-going muons K2K:2006odf ; MiniBooNE:2013qnd ;
MINOS:2014axb . One of the sources of these discrepancies was thought to be
coming from neutrino–nucleus interaction channels, which were not implemented
in the simulations used. The most-probable candidate is now believe to be the
so-called _multi-nucleon_ interaction, of which a similar process is known to
exist in electron–nucleus scattering. Inclusion of this interaction into
neutrino–nucleus simulations was discussed by Marteau in 1999 Marteau:1999kt .

In , the Valencia model by Nieves  _et al._ Gran:2013kda  is implemented. This
model considers an interaction involving the production of two nucleons and
two holes in a ground-state nuclear target (often called a 2p2h interaction in
the literature). Their model includes processes involving the exchange of
mesons between two nucleons and thus, sometimes, it is often referred to as a
Meson Exchange Current, or MEC, model. The model is not applicable for large
momentum transfer and thus the three-momentum transfer to the nucleus is
limited to $q_{3}<1.2$ GeV/c. This model does not predict how the four-
momentum is distributed between the two final-state nucleons, follows the
implementation in NuWro NuWro2006 . The directions of the outgoing nucleons
are selected to be uniformly distributed in the center of mass frame of the
nucleons. A separation energy (sometimes ambiguously referred to as the
_binding energy_) is subtracted from the energy transfer from the lepton
system. The outgoing nucleon momenta are required to be larger than the local
Fermi surface at the interaction position within the nucleus. The model for
the binding energy is described in Ref. Bourguille:2021:04  and depends on the
interaction position in the nucleus. For oxygen, typically between $\sim$50
MeV and $\sim$75 MeV of energy is lost to the nuclear response.

Neutral current multi-nucleon scattering is not yet implemented in .

####  3.2.3 Single Meson and Gamma Productions

Single pion production is one of the dominant neutrino interaction channels in
the few-GeV energy region. Therefore, it is important to understand these
interactions to study neutrino oscillations using the atmospheric or the
accelerator neutrinos. The particles from the single meson productions are
similar to the ones from nucleon decay, and thus, these interactions are also
important in the nucleon decay searches.

Single pion production in is implemented following the model by Rein and
Sehgal Rein:1980wg . An improved model, which takes into account the lepton
mass correction, by Berger and Sehgal Berger:2007rq  is also implemented. Both
of the models simulate these interactions in two steps. First, a neutrino
excites the nucleon and produce intermediate baryon resonance state, which
then decays into a single meson or gamma and baryon. The range of the
invariant mass of the intermediate state, $W$, is between the pion-production
threshold and smaller than 2 $\rm{GeV/c^{2}}$. In these two models, two kinds
of form factors are implemented. The first is a simple dipole form, as used in
the original publications by Rein and Sehgal Rein:1980wg . The axial vector
mass of the dipole type form factor was selected to be $1.21\rm{GeV/c^{2}}$
based on fits to K2K data. The second form was formalized by Graczyk and
Sobczyk Graczyk:2007xk  based on the Rarita-Schwinger formalization. The free
parameters of Graczyk-Sobczyk form were extracted by fits to hydrogen and
deuterium bubble chamber data. The direction of the pions in the resonance
rest frame, the so-called Adler frame, can be determined using the full
prescription by Rein Rein:1987cb  for all implemented resonances.
Alternatively, the Rein calculation can be only used for interactions
producing an intermediate $\Delta(1232)$ while, for higher-order resonances,
final-state pions are distributed uniformly in the Adler frame. For nuclear-
target interactions, Pauli-blocking is considered and final-state nucleon
produced in the resonance decay is required to have the momentum larger than
the local Fermi surface momentum modelled by an LFG. The overall effect is
small, typically less than a few percent of the events are rejected.

Single Kaon, Eta and $\gamma$ productions are simulated using the same
framework for the single pion production. The main differences are the decay
probabilities (branching ratios) of each simulated resonance to the relevant
final state. The production of other mesons, such as the omega, is not
simulated with the model described in this section. However, such particles
can be produced during the hadronization simulation in the Deep Inelastic
Scattering reactions, described in the next section.

For neutrino interactions with nucleon bound within a nucleus, produced
hadrons undergo the final state interactions as described in Sec. 3.3.
Additionally, formation zone effects are taken into account as described in
Sec. 3.2.6.

####  3.2.4 Shallow and Deep Inelastic Scattering

Shallow and deep inelastic scattering processes are separated as shown in
Table 3 to avoid double counting single- and multi-meson production channels.

Table 3: Shallow and deep inelastic scattering implementation in . $W$ is the invariant mass of the intermediate hadron system. | 1 meson | more than 1 meson
---|---|---
$W<2\rm{GeV/c^{2}}$ | (covered by single | custom multi-pion
| meson production ) | production model
$W>2\rm{GeV/c^{2}}$ | Pythia v5.72, included in CERNLIB 2005

Interactions producing more than one meson is simulated with a custom multi-
pion-production model. When $W$ is larger than 2 $\rm{GeV/c^{2}}$, all the
interactions, which produced at least one meson are simulated by PYTHIA v5.72
Sjostrand:1993yb , included in CERNLIB 2005\. For both cases, same parton
distribution functions (PDFs) are used. The default PDF is a modified version
of GRV98 Gluck:1998xa  that is based on the Bodek and Yang model Bodek:2005de
. The multi-pion production cross section is a function of the pion
multiplicity, which itself is a function of $W$, that models the probability
to produce more than one pion, with $W<2$ GeV/$c^{2}$ region. The enforcement
of multiple final-state pions avoids overlap with the resonance single meson
production channels. The mean multiplicity of charged hadrons is estimated
from the results of Fermilab 15-foot hydrogen bubble chamber experiment
MPIMULTP  and the Big European Bubble Chamber (BEBC) experiment MPIMULTP_BEBC
. Two parameter sets are available. The first uses the data from MPIMULTP  and
predicts

$\langle\textrm{n}_{\textrm{ch}}\rangle=0.09+1.83\times\mbox{ln}(W^{2}),$ (5)

for all channels. The second uses the fit results from both MPIMULTP  and
MPIMULTP_BEBC . This time, the parameter sets obtained by Bronner
Bronner:2016gmz  treat each combination of $\nu+p$, $\nu+n$, $\bar{\nu}+p$ and
$\bar{\nu}+n$ separately by fitting the measured average charged hadron
multiplicity data for each channel independently. The Bronner model for the
$\nu+p$ channel is reproduced in Fig. 6.

Figure 6: Results of the fit of the averaged charged hadron multiplicities as
a function of $W$ using the bubble chamber data sets MPIMULTP ; MPIMULTP_BEBC
. Figure reproduced from Ref. Bronner:2016gmz .

The obtained parameter sets are

$\displaystyle\langle\textrm{n}_{\textrm{ch}}\rangle$ $\displaystyle=$
$\displaystyle 0.58+1.35\times\mbox{ln}(W^{2})\leavevmode\nobreak\
\textrm{for}\leavevmode\nobreak\ \nu+p,$ (6)
$\displaystyle\langle\textrm{n}_{\textrm{ch}}\rangle$ $\displaystyle=$
$\displaystyle 0.35+1.24\times\mbox{ln}(W^{2})\leavevmode\nobreak\
\textrm{for}\leavevmode\nobreak\ \nu+n,$ (7)
$\displaystyle\langle\textrm{n}_{\textrm{ch}}\rangle$ $\displaystyle=$
$\displaystyle 0.41+1.18\times\mbox{ln}(W^{2})\leavevmode\nobreak\
\textrm{for}\leavevmode\nobreak\ \bar{\nu}+p,$ (8)
$\displaystyle\langle\textrm{n}_{\textrm{ch}}\rangle$ $\displaystyle=$
$\displaystyle 0.80+0.94\times\mbox{ln}(W^{2})\leavevmode\nobreak\
\textrm{for}\leavevmode\nobreak\ \bar{\nu}+n.$ (9)

For such interaction, we assume KNO scaling to determine the value of $W$
MPIWDIST . The forward-backward asymmetry of the charged hadron multiplicity
in the hadronic center of mass system is modelled as

$\frac{\textrm{n}_{\textrm{ch}}^{F}}{\textrm{n}_{\textrm{ch}}^{B}}=\frac{0.35+0.41\times\mbox{ln}(W^{2})}{0.5+0.09\times\mbox{ln}(W^{2})},$
(10)

which was derived from BEBC data MPIFBASYM .

The cross-section for the events with $W$ larger than 2 $\rm{GeV/c^{2}}$ is
calculated without the multiplicity factor because there is no overlap with
other implemented models. The kinematics of the produced particles are
determined by PYTHIA.

For neutrino interactions with a nucleon bound within a nucleus, final state
hadron interactions and formation zone effects are taken into account as for
single meson production.

####  3.2.5 Coherent and Diffractive Pion Productions

There are two coherent neutrino–nucleus pion production models implemented in
. The default model is based on the prescription by Berger and Sehgal
Berger:2008xs . This is an update to the previously publication model by Rein
and Sehgal Rein:1982pf , which is also implemented in . The new model uses
improved elastic pion–carbon cross-section data and lepton-mass effects are
properly taken into account. With these improvements, the model is applicable
to neutrinos with energies below several GeV.

A similar process, diffractive pion production, involves coherent-like pion
production, but with a single proton. implements the prescription by Rein
Rein:1986cd .

Pions produced through these channels are not affected by the final state re-
scattering, which is described in the following section 3.3.1.

####  3.2.6 Formation zone

The idea of a _formation zone_ is implemented and the production positions of
hadrons in nucleus are shifted from the initial neutrino interaction position.
The implemented model in is based on SKAT data Baranov:1985mb ; Ammosov:2001 .
The production points of the hadrons for those interactions are shifted using
the formation length ($L_{FZ}$), where

$L_{FZ}=p/\mu^{2},$ (11)

$p$ is the momentum of the hadron and $\mu=0.08$(GeV/c2). The actual size of
the shift is determined as $L_{FZ}\times(-\log({\rm rand}[0,1])$, where
rand[0,1] is a random number from 0 to 1. The distribution of secondary hadron
production position is shifted further from the center of the nucleus, the
produced position of hadrons shifts to the outer region of the nucleus and
thus reduce the probability of FSI.

###  3.3 Final-state Hadronic Re-scattering

simulates the interactions of pions, kaons, etas, omegas, protons and
neutrons, produced via neutrino interactions or nucleon decay, within the
nucleus. In order to simulate these hadron interactions, a custom semi-
classical intranuclear cascade (INC) model is used, as in most other
neutrino–nucleus simulations. In , a hadron produced in the nucleus is tracked
step by step from the production point until the particle escapes from the
nucleus. The size of each step is fixed at 0.2 fm. At each step, it is decided
whether the particle has interacted or not using the mean free paths for the
modelled interaction channels. A Woods-Saxon nucleon density function and the
local Fermi-gas model are used to determine the interaction positions and
kinematics of the initial and final states.

####  3.3.1 Meson interactions in nucleus

Among the modelled mesons, pion is the most important in the analyses of SK
and T2K. The mean free paths for the pion interactions in the Delta region are
calculated following the prescriptions by Salcedo  _et al._ Salcedo:1987md .
Their model takes into account the in-medium correction of $\Delta$ self
energy and uses the local Fermi-gas model. Therefore, the obtained mean free
paths are expressed as functions of nuclear density and pion momentum. This
model is applicable for the pions with momenta smaller than 500 MeV/c. The
mean free paths for pions with momenta larger 500 MeV/c was motivated by fits
to pion–nucleon scattering data and are expressed as functions of pion
momentum only.

In total six types of pion scattering channels are defined: low momentum pion
quasi-elastic scattering, low momentum pion charge exchange interaction, pion
absorption, high momentum pion quasi-elastic scattering, high momentum pion
charge exchange interaction, high momentum pion inelastic scattering (pion
production). Energy-independent normalization factors were defined for each of
these six channels and fit to pion–nucleus scattering data from various
experiments by Pinzon Guerra  _et al._ Elder:2019 , covering pion energies up
to 2GeV. Fig. 7 shows the comparison between the pion scattering data and the
simulated results from and the other simulation. The kinematics of particles
after the interaction are determined using the results of the phase shift
analysis Rowe:1978  with medium corrections as suggested by Seki  _et al._
Seki:1983 .

Figure 7: Comparisons of the $\pi^{+}$ absorption and charge exchange in
Carbon cross-sections data and the simulated results of three simulation
programs, GEANT4 (solid line), FLUKA (dashed line), and (dotted line). The
data are taken from duet:piscat ; jones:piscat ; ashery:pioncx ;
Bellotti1973_abs:piscat ; Bellotti1973_cx:piscat  and ABS (red) and CX
(black). This figure is taken from Ref. Elder:2019 .

The interactions of kaons, etas and omegas are treated similarly as high
momentum pions. The mean free paths are provided as functions of momentum for
each meson and the local Fermi-gas model is used to determine the kinematics
when an interaction occurs.

For kaons, (quasi-)elastic scattering and charge exchange interactions are
simulated. The mean free paths of these kaon interactions are extracted from
the results of the phase shift analysis of kaon scattering data Hyslop:1992cs
; Martin:1977me . These results are also used to determine the kinematics of
outgoing particles.

For etas, the $\eta N\rightarrow\eta N,\eta N\rightarrow\pi N^{\prime}$, and
$\eta N\rightarrow\pi\pi N^{\prime}$ processes are implemented. These
interactions are simulated assuming that an eta produces an excited nucleon
state, which then decays to give final state. The $N(1540),N(1650)$
intermediate baryon resonances are considered. The production cross-section of
baryon resonances are given by the Breit-Wigner formula:

$\sigma(k)=\frac{\pi}{k^{2}}\cdot\frac{\Gamma_{\eta
N}\Gamma_{X}}{(W-M_{N}^{*})^{2}+\Gamma^{2}_{tot}/4},$ (12)

where $W$ and $M_{N}^{*}$ are the invariant mass of the intermediate baryon
resonance and the mean mass of the resonance respectively, $\Gamma_{tot}$ is
the total width of $N^{*}$, $\Gamma_{\eta N}$ is the partial width of
$N^{*}\rightarrow\eta N$, $\Gamma_{X}$ is the partial width to the final state
$X$, where $X$ is the final state meson, $\eta$, $\pi$ or $\pi\pi$. The
produced particles are ejected isotropically in the intermediate resonance
rest frame.

For omegas, the $\omega N\rightarrow\pi N,\omega N\rightarrow\rho N,\omega
N\rightarrow\rho\pi N^{\prime},\omega N\rightarrow\rho\pi\pi N^{\prime},\omega
N\rightarrow\omega N,$ and $\omega N\rightarrow\sigma(f_{0})N^{\prime}$
processes are considered. The total and differential cross-sections of each
channel were calculated following the prescription by Lykasov  _et al._
Lykasov:1998ma .

####  3.3.2 Nucleon interactions in nucleus

The implementation of nucleon scattering is also based on the INC model. Three
types of interaction are considered, elastic scattering and one or two pion
production. These are implemented following the work by Bertini  _et al._
Bertini:1972vz , for MECC-7. The same nuclear density function and local
Fermi-gas model is used as the meson interactions in nucleus. When calculating
the interaction kinematics, an effective nucleon mass ($M_{N}^{eff}$) is used
instead of the free mass ($M_{N}^{free}$). Here, $M_{N}^{eff}$ is defined as

$M_{N}^{eff}=\sqrt{(M_{N}^{free}-8MeVc^{2})^{2}-(p_{F}^{surf})^{2}},$ (13)

where $P_{F}^{surf}$ is the density-dependent local Fermi surface momentum.
The momentum of the nucleon after an interaction is required to be larger than
the local $P_{F}^{surf}$, which is the standard Pauli-blocking procedure in
the local Fermi-gas model framework. Produced pions are assumed to arise from
the decay of $\Delta(1232)$ resonances produced by nucleon- nucleon
scattering. To simulate the $\Delta(1232)$ production a simple isobar model
Lindenbaum:1957ec  is used. When a pion is produced in the nucleus after the
nucleon re-scattering, that pion is independently tracked from the point of
generation using the pion transport simulation described in section 3.3.1. The
predicted cross-sections for proton-carbon scattering is shown in Fig. 8.

Figure 8: Proton-carbon scattering cross-sections as a function of the proton
kinetic energy, as predicted by . To obtain this plot, protons are injected
from outside of the nucleus and simulates the interactions in the carbon
nucleus. The black line shows the total cross-section, the purple line shows
the elastic scattering cross-section and the blue line shows the inelastic
scattering cross-sections.

###  3.4 Known Implementation Limitations

is primarily developed to enable world-leading SK and T2K neutrino
oscillation, proton decay, and neutrino cross-section measurements, for which
a fully consistent description of Nature is desirable but not necessary, and
as such has some known limitations beyond the imperfect predictions of the
implemented physics models.

The most prominent limitation is that the modelling of the initial nuclear
state is tied to the individual interaction channel being simulated. For most
channels, the initial state is always modelled as a Fermi-gas, for multi-
nucleon scattering an LFG is used, and only for quasi-elastic scattering is
the more sophisticated spectral function an option. The SF is only designed to
model QE interactions, and only for specific nuclei where relevant electron-
scattering data exist, for other nuclear targets a Fermi-gas is used. While,
as noted above, the modelling of hadron re-scattering is based on the density
and momentum predictions from an LFG model. Such an inconsistent model is
sometimes affectionately referred to as a _Franken-model_ , after the
fictional scientist and his Gothic horror implementation. For single meson
production, nuclear effects beyond the initial-state nuclear momentum and
Pauli-blocking are ignored. It might be expected that the inclusion of the
initial-state nucleon removal energy would affect the predictions. As
statistical uncertainties are reduced with the next generation of long-
baseline oscillation experiments, addressing these inconsistencies will become
a focus of future development.

does not attempt to model a number of interaction channels that are irrelevant
to SK and T2K analyses, these include neutrino–electron elastic scattering
(used to provide complementary neutrino flux normalization constraints for
higher energy beams PhysRevD.100.092001 ) and inverse beta decay (important
for simulating reactor neutrino experiments).

##  4 Additional Tools

So far we have focused on the implementation and physics models for simulation
single neutrino–nucleus interactions. While this constitutes the core of the
simulation, to be an effective tool for data analysis, provides a number of
other tools and features.

#####  neutgeom:

The neutgeom tool is used to interface to a neutrino beamline simulation,
perform neutrino ray-tracing, interface with a detector geometry description
and correctly distribute interaction positions through the detector.
Currently, neutgeom is only able to consume neutrino flux _vectors_ from
JNUBEAM, which simulation the J-PARC neutrino beam for the T2K experiment.
There is some interest in distributing neutgeom as a standalone tool capable
of interfacing with other neutrino beam simulations and interaction
simulations.

#####  NReWeight:

Cross section _reweighting_ is an important tool for systematic error
estimation for neutrino-scattering analyses. For systematic parameters that
can be effectively reweighted it enables interaction model variations to be
applied at analysis time, rather than requiring re-simulation and re-analysis,
reducing the computing time taken to estimate systematic uncertainties by many
many orders of magnitude. At its core, the reweighting technique calculates a
weight for each already-simulated interaction,

$W=\frac{\sigma\left(p^{\prime},\vec{x}\right)}{\sigma\left(p,\vec{x}\right)},$

where $\vec{x}$ encapsulates the kinematics of the simulated interaction and
$p$ and $p^{\prime}$ fully describe the model choices and any free parameters
in the already-simulated and the varied model, respectively. This procedure is
only exact when the varied model, $p^{\prime}$, predicts a range of $\vec{x}$
that is the same or a strict subset of the range predicted by $p$. provides
the NReWeight package that implements exact cross section reweighting enabling
the variation of nucleon form factors for CCQE and single meson production
parameters after simulation. NReWeight also exposes reweighting for the meson
and nucleon intranuclear cascade, which is exact for modest variations of the
underlying meson–nucleon and nucleon–nucleon scattering probabilities.
NReWeight is a critical tool for SK and T2K analyses.

##### Interface to GEANT3 and GEANT4:

The simulation of an intranuclear cascade and of meson–nuclear and
nucleon–nuclear scattering is conceptually similar, the difference is whether
the probe beings the simulation inside or outside of the simulated nucleus. As
described in section 3.3.1, the pion–nucleon interaction cross-sections were
tuned to pion–nuclear scattering data. As a result, it is attractive for
physics consistency to be able to simulate pion–nucleus scattering with the
INC model. To achieve this, an interface between and GEANT was developed and
integrated to the T2K near and SK detector simulation programs.

##  5 Future direction

Current and future neutrino experiments require precise neutrino interaction
simulations to achieve their ambitious physics goals. To meet these
requirements, will be continue to be developed and improved.

Current modelling improvements include the implementation of the state-of-the-
art CCQE and multi-nucleon model by Amaro  _et al._ Amaro:2021sec , and single
pion production models by Kabirnezhad Kabirnezhad:2020wtp  and by Sato  _et
al._ Sobczyk:2018ghy .The implementation of a rudimentary QE-only electron-
scattering simulation in is underway. This will enable improved extended
validations of the implemented physics that is most critical to T2K analyses.

The dependence on CERNLIB is problematic, as the library is no longer
maintained. Building and distributing the library for modern compilers and
operating systems takes time away from more important development efforts.
currently relies on CERNLIB for reading configuration files, random number
generation and some common mathematical operations, and PYTHIA v5.72 for
simulating SIS and DIS interactions. Resolving the dependence on PYTHIA v5 is
not simple as changing to a newer version (v6 or v8) would require the
implementation of the particles’ kinematics determination (vector generation)
functions in , as the PYTHIA implementation that we rely on has been dropped
or does not cover the entire kinematic region. Removing the dependency on
CERNLIB is a high priority for near-future maintenance.

The current closed source nature of is undesirable. Exposure to more users and
use cases will result in code, interface, and physics improvements. However,
the lack of human resources render it difficult to support as a more general
tool. Work has begun, in collaboration with other neutrino interaction
simulation stake-holders, to define, test, and implement a new community-
designed event format and event generation API barrow2020summary . These
critical future developments will be implemented in as they become defined and
mature.

##  6 Summary

is a general purpose neutrino interaction simulation used and improved by
members of the SK and the T2K collaborations, with critical additional
contributions from interaction theory groups. Future development will target
the physics requirements of Super-Kamiokande, T2K, and Hyper-Kamiokande, and
software integration improvements to the API and data formats.

##  7 Acknowledgments

The authors would like to thank the many contributors of , and the users in
the SK, T2K, Ninja, and WAGASCI collaborations that push us to continue to
develop it as a world-class analysis tool. The authors would also like to
thank Kendall Mahn for help proofreading this document. This work is partially
supported by MEXT/JSPS KAKENHI Grant Number 18H05536. This material is based
upon work supported by the U.S. Department of Energy, Office of Science,
Office under Award Number(s) DE-SC0015903.

## References

  * (1) E. S. Pinzon Guerra et al., Phys. Rev. D 99, 052007 (2019)
  * (2) M. Nakahata et al. [Kamiokande], J. Phys. Soc. Jap. 55, 3786 (1986)
  * (3) Y. Fukuda et al. [Super-Kamiokande], Phys. Rev. Lett. 81, 1562 (1998)
  * (4) K. Abe et al. [Super-Kamiokande], Phys. Rev. D 97, no.7, 072001 (2018)
  * (5) K. Abe et al. [T2K], Phys. Rev. D 91, no.7, 072010 (2015)
  * (6) K. Abe et al. [T2K], Phys. Rev. D 87, no.1, 012001 (2013)
  * (7) M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara and S. Midorikawa, Phys. Rev. D 92, no.2, 023004 (2015)
  * (8) O. Benhar, and N. Farina, and H. Nakamura, and M. Sakuda, and R. Seki, Phys. Rev. D 72, 053005 (2005)
  * (9) H. Ejiri, Phys. Rev. C 48, 1442 (1993)
  * (10) S. Jadach, Z. Was, R. Decker and J. H. Kuhn, Comput. Phys. Commun. 76, 361-380 (1993)
  * (11) C.H.Llewellyn Smith, Phys. Rep. 3C, 261 (1972)
  * (12) R.A. Smith and E.J. Moniz, Nucl. Phys. B 43, 605 (1972)
  * (13) R. Gran, J. Nieves, F. Sanchez and M. J. Vincente Vacas, Phys. Rev. D 88, 113007 (2013)
  * (14) B. Bourguille, and J. Nieves, and F. Sanchez, JHEP 04, 004 (2021)
  * (15) J. Nieves, J. E. Amaro and M. Valverde, Phys. Rev. C 70, 055503 (2004) [erratum: Phys. Rev. C 72, 019902 (2005)]
  * (16) O. Benhar, A. Fabrocini, S. Fantoin and I. Sick, Nucl. Phys. A 579, 493 (1994)
  * (17) C. Juszczzak, J. A. Nowak, and J. T. Sobczyk, Nucl. Phys. B Proceedings Supplements, 159, 211 (2006)
  * (18) A. Furmanski, University of Warwick (2015), PhD Thesis No. 74036
  * (19) R. Bradford, A. Bodek, H. S. Budd and J. Arrington, Nucl. Phys. B Proc. Suppl. 159, 127 (2005)
  * (20) A. Bodek, S. Avvakumov, R. Bradford and H. S. Budd, Eur. Phys. J. C 53, 349 (2008)
  * (21) B. Bhattacharya, R. J. Hill and G. Paz, Phys. Rev. D 84, (2011), 073006.
  * (22) A. S. Meyer, and M. Betancourt, and R. Gran, and R. J. Hill, Phys. Rev. D 93, (2016), 113015.
  * (23) C. Adamuscin, E. Tomasi-Gustafsson, E. Santopinto and R. Bijker, Phys. Rev. C 78, 035201 (2008)
  * (24) J. Stowell, University of Sheffield (2019), PhD Thesis No. 22623
  * (25) R. Gran et al. [K2K], Phys. Rev. D 74, 052002 (2006)
  * (26) A. A. Aguilar-Arevalo et al. [MiniBooNE], Phys. Rev. D 88, no.3, 032001 (2013)
  * (27) P. Adamson et al. [MINOS], Phys. Rev. D 91, no.1, 012005 (2015)
  * (28) J. Marteau, Eur. Phys. J. A 5, 183 (1999)
  * (29) D. Rein and L. M. Sehgal, Annals Phys. 133, 79 (1981)
  * (30) C. Berger and L. M. Sehgal, Phys. Rev. D 76, 113004 (2007)
  * (31) K. M. Graczyk and J. T. Sobczyk, Phys. Rev. D 77, 053003 (2008)
  * (32) D. Rein, Z. C Phys. 35, 43 (1987)
  * (33) T. Sjostrand, Comput. Phys. Commun. 82, 74- (1994)
  * (34) M. Glück, E. Reya and A. Vogt, Eur. Phys. J. C 5, 461-470 (1998)
  * (35) A. Bodek and U. k. Yang, AIP Conf. Proc. 792, no.1, 257 (2005)
  * (36) S.J. Barish et al., Phys. Rev. D 17, 1 (1978)
  * (37) B. JongeJans et al., Nuovo Cim. A 101, 435 (1989)
  * (38) C. Bronner and M. Hartz, JPS Conf. Proc. 12, 010041. (2016)
  * (39) H. Saarikko, Proceedings of 9th International Neutrino Conference, 507 (1979)
  * (40) S. Barlag et al. Z.Phys. C 11, 283 (1982)
  * (41) C. Berger and L. M. Sehgal, Phys. Rev. D 79, 053003 (2009)
  * (42) D. Rein and L. M. Sehgal, Nucl. Phys. B 223, 29 (1983)
  * (43) D. Rein, Nucl. Phys. B 278, 61 (1986)
  * (44) D. S. Baranov et al.,Yad. Fiz.,40, 1454 (1984)
  * (45) V. Ammosov, Low multiplicity final states in neutrino interaction at low energy – SKAT results, NuINT01 workshop, KEK, Tsukuba, Japan (2001)
  * (46) L. L. Salcedo, E. Oset, M. J. Vicente-Vacas and C. Garcia-Recio, Nucl. Phys. A 484, 557-592 (1988)
  * (47) E. S. Pinzon Guerra et al. (DUET Collaboration), Phys. Rev. C 95, 045203. (2017)
  * (48) M. K. Jones et al., Phys. Rev. C 48, 2800 (1993)
  * (49) D. Ashery et al., Phys. Rev. C 30, 946 (1984)
  * (50) E. Bellotti, D. Cavalli, and C. Matteuzzi, Nuovo Cim. A 18, 75 (1975)
  * (51) E. Bellotti, S. Bonetti, D. Cavalli, and C. Matteuzzi, Nuovo. Cim. A 14 (1973), 567.
  * (52) G. Rowe, M. Salomon, and R. H. Landau, Phys. Rev. C 18, 584 (1978)
  * (53) R. Seki, M. Masutani, Phys. Rev. C 27, 2799 (1983)
  * (54) J. S. Hyslop, R. A. Arndt, L. D. Roper and R. L. Workman, Phys. Rev. D 46, 961 (1992)
  * (55) B. R. Martin and M. K. Pidckco Nucl. Phys. B 126, 285 (1977)
  * (56) G. I. Lykasov, W. Cassing, A. Sibirtsev and M. V. Rzyanin, Eur. Phys. J. A 6, 71 (1999)
  * (57) H. W. Bertini, Phys. Rev. C 6, 631 (1972)
  * (58) S. J. Lindenbaum and R. M. Sternheimer, Phys. Rev. 105, 1874 (1957)
  * (59) Valencia, E. et al. [MINER$\nu$A], Phys. Rev. D 100, no.9, 092001 (2019)
  * (60) J. E. Amaro, M. B. Barbaro, J. A. Caballero, T. W. Donnelly, R. Gonzalez-Jimenez, G. D. Megias and I. R. Simo, [arXiv:2106.02857 [hep-ph]]
  * (61) M. Kabirnezhad, Phys. Rev. D 102, no.5, 053009 (2020)
  * (62) J. E. Sobczyk, E. Hernández, S. X. Nakamura, J. Nieves and T. Sato, Phys. Rev. D 98, no.7, 073001 (2018)
  * (63) J. Barrow  _et al._ , arxiv [hep-ex] 2008.06566, FERMILAB-CONF-20-423-ND-SCD