\subsection{Changing the $\pi^+$ inelastic cross section}
Figure \ref{fig:pion_inel} shows comparison plots for the following samples: 1) running Geant4 with both of the nominal cross sections (``Nominal'' in the plots), 2) running Geant4 with the inelastic cross section multiplied by 1.5 (``Varied''), and 3) applying Geant4Reweight to the nominal cross section sample according to the inelastic cross section change (``Weighted''). In figure \ref{subfig:pion_inel_len}, the track length of the particle is shown. This is the integrated path length traveled by the $\pi^+$ from its creation to the point at which it stops, decays in flight, or undergoes an inelastic interaction. Increasing the inelastic cross section increases the chance that it will interact earlier, before stopping or decaying, thus shifting the Weighted and Varied distributions lower, as expected. Figure \ref{subfig:pion_inel_nElast} shows the number of elastic scatters performed by the $\pi^+$ during its travel. Similarly, increasing the inelastic cross sections means the $\pi^+$ will be less likely to travel as far as in the nominal simulation, suppressing the chance for elastic scatters to occur. The figures show good agreement between the varied and weighted samples, supporting the validity of our reweighting scheme.

\begin{figure}[htbp]
    \centering
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/len_1_14_21_bias.pdf}
        \caption{}
        \label{subfig:pion_inel_len}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.45\textwidth}
        \centering

        \includegraphics[width=\textwidth]{Figures/nElast_1_14_21_bias.pdf}
        \caption{}
        \label{subfig:pion_inel_nElast}
    \end{subfigure}
    \caption{Comparison of pion track lengths (left) and number of elastic scatters (right) using Geant4 with nominal cross sections (Nominal) and with inelastic cross sections varied directly in Geant4 (Varied) and using Geant4Reweight (Weighted).  }
    \label{fig:pion_inel}
\end{figure}

\subsection{Changing the $\pi^+$ elastic cross section}
Figure \ref{fig:pion_el} shows similar comparison plots as the previous section, but for the sample with the elastic cross section multiplied by 1.5 and the inelastic cross section kept the same. In figure \ref{subfig:pion_el_len}, the track length distribution shows no significant change. This is due to the track-length definition above, as well as the fact that the elastic interaction in Geant4 removes minimal energy from the $\pi^+$. As such, more elastic scatters do not cause the $\pi^+$ to stop earlier. Figure \ref{subfig:pion_el_nElast} shows that the number of elastic scatters performed by the $\pi^+$ increases as expected.

\begin{figure}[htbp]
    \centering
    \begin{subfigure}[t]{0.45\textwidth}
        \centering

        \includegraphics[width=\textwidth]{Figures/len_1_14_21_elast.pdf}
        \caption{}
        \label{subfig:pion_el_len}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.45\textwidth}
        \centering

        \includegraphics[width=\textwidth]{Figures/nElast_1_14_21_elast.pdf}
        \caption{}
        \label{subfig:pion_el_nElast}
    \end{subfigure}
    \caption{Comparison of pion track lengths (left) and number of elastic scatters (right) using Geant4 with nominal cross sections (Nominal) and with elastic cross sections varied directly in Geant4 (Varied) and using Geant4Reweight (Weighted).  }
    \label{fig:pion_el}
\end{figure}

\subsection{Recovering the nominal cross section}
As a final demonstration, we show the ability for Geant4Reweight to recover the nominal cross sections in Geant4 by reweighting the samples generated with varied cross sections. This is shown first in figure~\ref{fig:pion_inel_2}. Here, ``Nominal'' and ``Varied'' are samples generated with nominal cross sections and with the inelastic cross section scaled by 1.5. ``Weighted*'' shows a reweighting applied to the Varied sample, where a scaling factor of 2/3 is applied to the varied cross section to return to the nominal. The figure shows the ability to reweight from a varied cross section back to nominal.

Figure \ref{fig:pion_el_2} shows a similar treatment, but with the elastic cross section increased by 50\% and then weighted to recover the nominal model. The nominal elastic cross section can be recovered via reweighting as shown.

\begin{figure}[htbp]
    \centering
    \begin{subfigure}[t]{0.45\textwidth}
        \centering

        \includegraphics[width=\textwidth]{Figures/len_1_14_21_revert.pdf}
        \caption{}
        \label{subfig:pion_inel_len_2}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.45\textwidth}
        \centering

        \includegraphics[width=\textwidth]{Figures/nElast_1_14_21_revert.pdf}
        \caption{}
        \label{subfig:pion_inel_nElast_2}
    \end{subfigure}
    \caption{Comparison of pion track lengths (left) and number of elastic scatters (right) using Geant4 with nominal cross sections (Nominal) and with inelastic cross sections varied directly in Geant4 (Varied) and using Geant4Reweight applied to the Varied sample  (Weighted*).}
    \label{fig:pion_inel_2}
\end{figure}

\begin{figure}[htbp]
    \centering
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/len_1_14_21_revert_elast.pdf}
        \caption{}
        \label{subfig:pion_el_len_2}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/nElast_1_14_21_revert_elast.pdf}
        \caption{}
        \label{subfig:pion_el_nElast_2}
    \end{subfigure}
    \caption{Comparison of pion track lengths (left) and number of elastic scatters (right) using Geant4 with nominal cross sections (Nominal) and with elastic cross sections varied directly in Geant4 (Varied) and using Geant4Reweight applied to the Varied sample  (Weighted*).}
    \label{fig:pion_el_2}
\end{figure}

\subsection{Parameter definitions}
The parameters used within the fit are user-defined variations across regions of incident hadron momentum, for the various cross sections. For charged pions, this includes the exclusive cross sections defined in table~\ref{tab:pi_exc}, as well as the total inelastic (the sum of the exclusive channels --- also referred to as reaction) and elastic cross sections. The reaction cross section parameters work simultaneously with the exclusive parameters, similar to the weighting scheme in section~\ref{sec:reweighting}. Flat variations are applied across the momentum range chosen by users.

During each iteration of the fit, the values of the user-defined parameters are changed. The cross sections are scaled according to the variations defined by the new parameter values, and fit to the provided data. If the fit is successful, the best-fit values of the parameters and their covariances are saved in the output ROOT file.

\subsection{$\chi^2$ definitions}

The data in the fit is categorized by the interaction channel of each measurement. The relevant channels are the exclusive channels (such as those listed in table \ref{tab:pi_exc}), combinations of these (such as pion absorption + charge exchange), the elastic interaction, the reaction (total inelastic) interaction, and the inclusive interaction (reaction + elastic). The quantity of measurements for each exclusive channel varies widely, and a simple $\chi^2$ calculation will likely lead to channels with more measurements dominating the fit. Instead, a reduced $\chi^2$ statistic is calculated for each channel, as shown in equation \ref{eq:exc_chi2}. Here, $k$ denotes the interaction channel, $N_k$ is the number of cross section measurements for channel $k$, and $p_k$ is the number of user-defined parameters for that channel. $\sigma^{\textrm{MC}}$ denotes the Geant4 cross section that has been varied by the fit parameters as described above.

\begin{align}
\chi^2_k &= \frac{1}{N_k - p_k} \sum_{i}^{N_k} \left(\frac{\sigma^{\textrm{
Data}}_i - \sigma^{\textrm{MC}}_i}{\Delta\sigma^{\textrm{Data}}_i} \right)^2\label{eq:exc_chi2} \\
\chi^2 &= \frac{1}{N_{\textrm{C}} + f - 1}\bigg(\sum\limits_{k} \chi^2_k + f\chi^{2}_{\textrm{Total}}\bigg)  \label{eq:total_chi2}
\end{align}
The full $\chi^2$ is shown in equation \ref{eq:total_chi2}. Here, the $\chi^2$ contributions from all channels other than the inclusive interaction are simply summed. The inclusive interaction is treated slightly differently because variations to the total cross section can sometimes dominate the fit. As such, a scaling parameter $f$ is introduced, allowing the user to control the contribution of the inclusive channel to the full $\chi^2$. The leading factor is the effective number of degrees of freedom, where $N_C$ is the number of channels included in the sum. This term ensures $\chi^2 \sim 1$ as $\big(\sigma^{\textrm{
Data}}_i - \sigma^{\textrm{MC}}_i\big)$ approaches  $\Delta\sigma^{\textrm{Data}}$.

\subsubsection{Fit with correlated data}
The reduced $\chi^2$ statistic defined in equation \ref{eq:exc_chi2} is valid for uncorrelated data.
However, a recent measurement by the DUET collaboration includes correlations in their published results for $\pi^+$ absorption and charge exchange on carbon~\cite{DUET}. We have included this data in the fit in a specific way to account for these correlations. A standard $\chi^2$ definition is shown in equation \ref{eq:DUET_chi2}. For this set of data, the absorption and charge exchange cross sections were each measured at five momenta. As such, the sum in this equation runs over all ten data points. In order to treat the measurements of these two channels in a similar manner as above, we have chosen to split this value equally between the modified reduced $\chi^2$ statistics for the absorption and charge exchange channels, as shown in equations \ref{eq:abs_chi2_special} and \ref{eq:cex_chi2_special}. The denominator of the first term in each of these equations includes the five data points for each channel.

\begin{equation}\label{eq:DUET_chi2}
    \chi^2_{\textrm{DUET}} = \sum\limits_{i,j}^{10}\big(\sigma^{\textrm{DUET}}_i - \sigma^{\textrm{MC}}_i\big)V_{ij}^{-1} \big(\sigma^{\textrm{DUET}}_j - \sigma^{\textrm{MC}}_j\big)
\end{equation}

\begin{equation}\label{eq:abs_chi2_special}
    \chi^2_{\textrm{Abs}} = \frac{1}{N_{\textrm{Abs}} + 5 - p_{\textrm{Abs}}}\Bigg(\sum\limits^{N_{\textrm{Abs}}}_{i}\bigg(\frac{\sigma_i^{\textrm{Data}} - \sigma_i^{\textrm{MC}}}{\Delta\sigma^{\textrm{Data}}_i}\bigg)^2 + \frac{\chi^2_{\textrm{DUET}}}{2}\Bigg)
\end{equation}

\begin{equation}\label{eq:cex_chi2_special}
    \chi^2_{\textrm{Cex}} = \frac{1}{N_{\textrm{Cex}} + 5 - p_{\textrm{Cex}}}\Bigg(\sum\limits^{N_{\textrm{Cex}}}_{i}\bigg(\frac{\sigma_i^{\textrm{Data}} - \sigma_i^{\textrm{MC}}}{\Delta\sigma^{\textrm{Data}}_i}\bigg)^2 + \frac{\chi^2_{\textrm{DUET}}}{2}\Bigg)
\end{equation}
In general, other correlated data can be included in the fit in a similar manner. If such a set of data exists, a user can implement their own calculation of its contribution to the specific $\chi^2_k$.

\subsection{Example fit}For demonstration, we have fit Geant4 to a collection of $\pi^+$-C data. The parameters we have chosen to use are shown in table \ref{tab:parameters}. The data included in the fit is shown in table \ref{tab:fit_data}. Throughout this section, the ``reaction'' channel is used as another term for any inelastic pion-nucleus interaction. Its cross section is the combination of all exclusive channels (i.e. it is the total inelastic cross section shown in figure \ref{fig:xsec}).

\begin{table}[htbp]
\centering
\caption{\label{tab:parameters} The parameters used within the demonstration fit.}
\smallskip
    \begin{tabular}{|l|l|l|}
    \hline
    Parameter Name & Channel & Momentum Range [MeV$/c$]\\
    \hline
    fReacLow & Reaction & 10--200 \\

    fReacHigh & Reaction & 700--2005 \\

    fAbs & Absorption & 200--700 \\

    fCex & Charge Exchange & 200--700 \\

    fInel & Quasielastic & 200--700 \\
    \hline
    \end{tabular}
\end{table}

\begin{table}[htbp]
\centering
\caption{\label{tab:fit_data} The data sets used within the demonstration fit.}
\smallskip
    \begin{tabular}{|l|l|l|}
    \hline
    Reference & Channel(s) & Momentum Range [MeV/$c$]\\
    \hline
    B. W. Allardyce \textit{et al.}~\cite{Allardyce} & Reaction & 710--2000\\

    D. Ashery \textit{et al.}~\cite{Ashery2173} & Quasielastic, Abs. + Ch. Ex. & 175--433\\

    D. Ashery \textit{et al.}~\cite{Ashery946} & Charge Exchange & 265\\

    E. Bellotti \textit{et al.}~\cite{Bellotti} & Absorption, Charge Exchange & 230\\

    S. M. Levenson \textit{et al.}~\cite{Levenson} & Quasielastic & 194--417\\

    O. Meirav \textit{et al.}~\cite{Meirav} & Reaction & 128--170\\

    I. Navon \textit{et al.}~\cite{Navon} & Abs. + Ch. Ex. & 128\\

    E. S. Pinzon-Guerra \textit{et al.}~\cite{DUET} & Absorption, Charge Exchange & 206--295 \\

    A. Saunders \textit{et al.}~\cite{Saunders} & Reaction & 116--150\\
    \hline
    \end{tabular}
\end{table}

The results are shown in figures \ref{fig:pars_cov}--\ref{fig:scans}. Figure \ref{fig:pars_cov} shows the best fit values for the parameters and their covariance. Of note are the slightly larger covariances between the absorption and charge exchange parameters which are introduced due to including the combined absorption and charge exchange cross section data. In figure \ref{fig:fit}, cross section data used in the fit is compared to the nominal Geant4 prediction (blue line) and the best fit (black line) \& $\pm1\sigma$ error bands (red area). These figures highlight a current limitation in the framework as seen in the discontinuities in the best fit values of the cross sections and the error band sizes at the parameter transitions. Highly precise reaction data is included for momentum above 700 MeV/$c$. As such, this reduces the error on the high-momentum parameter as compared to the lower-momentum regions on the parameter affecting that region when compared to the middle region, which is fit to less accurate data.
Finally, in figure \ref{fig:scans} are 1-D parameter scans of the fit $\chi^2$ around the best-fit point. All plots are generated by the fitting application and included in the output.

\begin{figure}[htbp]
    \centering
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/parameters.pdf}
        \caption{Best-fit parameter values and 1-D uncertainty.}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/covariance.pdf}
        \caption{Output covariance matrix.}
    \end{subfigure}
    \caption{}
    \label{fig:pars_cov}
\end{figure}

\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[t]{0.4\textwidth}
        \centering

        \includegraphics[width=\textwidth]{Figures/Fit/reaction.pdf}
        \caption{}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/inelastic.pdf}
        \caption{}
    \end{subfigure}

        \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/abs.pdf}
        \caption{}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/cex.pdf}
        \caption{}
    \end{subfigure}

    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/abscx.pdf}
        \caption{}
    \end{subfigure}
    \caption{Fit results of $\pi^+$-C cross sections extracted from Geant4 (with the Bertini Cascade model) fit to data (listed in table \ref{tab:fit_data}). The parameters used in the fit are defined in table \ref{tab:parameters}. The blue dashed line is the nominal cross section when extracted from Geant4, the black dashed line is the best fit value, and the red bands are calculated by setting the parameters to their $\pm1\sigma$ values.}
    \label{fig:fit}
\end{figure}

\begin{figure}[!htbp]
    \centering
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/fReacLow.pdf}
        \caption{fReacLow}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/fReacHigh.pdf}
        \caption{fReacHigh}
    \end{subfigure}

    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/fAbs.pdf}
        \caption{fAbs}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/fCex.pdf}
        \caption{fCex}
    \end{subfigure}

    \begin{subfigure}[t]{0.4\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/Fit/fInel.pdf}
        \caption{fInel}
    \end{subfigure}
    \caption{1-D parameter scans of the fit $\chi^2$ near the best-fit point.}
    \label{fig:scans}

\end{figure}

\clearpage

\section{Introduction}
\label{sec:intro}
\input{introduction}

\section{Geant4 simulation strategy}
\label{sec:geant4}
\input{transport}

\section{Reweighting technique}
\label{sec:reweighting}
\input{reweighting}

\section{Reweighting examples}
\label{sec:example}
\input{example}

\section{Fitting framework}
\label{sec:fitting}
\input{fitting}

\section{Summary}
Geant4Reweight will serve the neutrino physics community by providing a framework that allows users to estimate and propagate systematic uncertainties related to hadron--nucleus interactions in Geant4. This paper describes the framework and its validation. Multiple experiments --- including MicroBooNE, NOvA, and ProtoDUNE-SP --- have begun to use this framework for uncertainty propagation in their analyses.

\acknowledgments
This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institute for Science and Education for the DOE under contract number DE‐SC0014664. Authors are also supported by DOE Award Number DE-SC0015903. We acknowledge the hard work of the Geant4 collaboration in the creation, documentation, and maintenance of the Geant4 software package, without which this work would not be possible. We also gratefully acknowledge a careful reading of this manuscript and useful comments by Wayne Repko.

\subsection{Basis for weights}
The previous section provided a description of the Geant4 simulation strategy and introduced relevant concepts such as processes, cross sections, interaction lengths, and models in the parlance of Geant4's simulation. This section describes the reweighting technique used to create weights for particle tracks under some given change to a process' cross section. The processes we are interested in reweighting (the hadronic elastic and inelastic interactions) are \textit{discrete} processes, so we will focus the remaining discussion on this type of process. To inform the description of the reweighting process, it is helpful to construct event probabilities that can be varied according to the cross sections. The following discussion describes this.

When the particle track is first created, each \textit{discrete} process samples the exponential distribution:
\begin{equation}
p(N) = e^{-N}dN.
\end{equation}
This provides the process with a number of interaction lengths ($N$). The interaction length proposed for a step is $N\lambda$, where $\lambda$ is the mean free path of the process for that step and is inversely proportional to the cross section for that process. At each step, if the process in question did not occur, the process subtracts an amount equal to $n = \Delta L/\lambda$ from $N$ where $\Delta L$ is the step length taken and $n$ is the number of interaction lengths taken by the step. Once this process is invoked in the procedure detailed in section \ref{sec:geant4}, $N$ is sampled again. Thus, for a given process, the probability for that process to \textit{not} occur after a series of steps ($P_S$ in short for $P_{Survive}$) is given by the probability to sample $N$ greater than the sum of all interaction lengths traveled by particle ($N_T$):

\begin{equation}\
\begin{split}
   P_{S} &= P(N > N_T)
   \\
   &= \int_{N_{T}}^{\infty}e^{-N}dN = e^{-N_{T}} = e^{-\sum\limits_i^{f} n_{i}}
   \\
   &= e^{-\sum\limits_i^{f} \Delta L_i/\lambda_i} = e^{-\sum\limits_i^{f} \Delta L_i\sigma_i}.
\end{split}
\end{equation}
where $n_i = \Delta L_i/\lambda_i$ is the number of interaction lengths traveled in step $i$. The last term uses the cross section $\sigma_i = 1/\lambda_i$ to simplify notation.
Normally, $1/\lambda = \sigma\rho N_A/M$, where $M$ and $\rho$ are the molar mass and density of the surounding material, and $N_A$ is Avogadro's number. The extra factors have been absorbed into $\sigma$ for clarity.
Note that the cross section $\sigma_i$ is generally dependent on the particle's momentum and the surrounding material for step $i$. This probability can be rewritten as the product of the probabilities for each step to occur without interaction:
\begin{equation}\label{eqn:p_survive}
    P_{S} = \prod\limits_i^{f}P_{S, i} = \prod\limits_i^{f}e^{-\Delta L_i\sigma_i}.
\end{equation}
If the process does occur, it only occurred in the last step, so the probability ($P_I$ in short for $P_{Interact}$) is largely unchanged except for the term contributing to the final step:
\begin{equation}\label{eqn:p_interact}
    \begin{split}
    P_{I} &= \bigg(\prod\limits_i^{f - 1}P_{S,i}\bigg)\bigg(P_{I, f}\bigg) = \bigg(\prod\limits_i^{f - 1}P_{S,i}\bigg)\bigg(1 - P_{S, f}\bigg)
    \\
    &= \bigg(\prod\limits_i^{f - 1}e^{-\Delta L_i\sigma_i}\bigg)\bigg(1 - e^{-\Delta L_f\sigma_f}\bigg).
    \end{split}
\end{equation}

These probabilities can be easily extended to consider multiple processes. In equations \ref{eqn:p_survive} and \ref{eqn:p_interact}, the cross sections will be extended to be sums of the cross sections of all processes.
\begin{equation}\label{eq:total_sigma}
    \sigma = \sigma_I + \sigma_E + \sigma_X
\end{equation}
For our purposes, we are interested in reweighting the elastic and inelastic hadronic interactions, so the cross sections for these processes are explicitly stated in equation \ref{eq:total_sigma} as $\sigma_I$ and $\sigma_E$ respectively. All other processes are combined into a single effective cross section $\sigma_X$. The processes included in $\sigma_X$ are Coulomb scattering and decay (for the relevant particles) which is effectively treated as an interaction with an interaction length given by $1/(\gamma \beta c \tau)$. Thus, for a given process $p$, the probability it will occur is given by
\begin{equation}\label{eqn:p_interact_spec}
    P_{I, p} = \bigg(\prod\limits_i^{N_f - 1}e^{-\Delta L_i\sigma_i}\bigg)\bigg(1 - e^{-\Delta L_f\sigma_f}\bigg)\bigg(\frac{\sigma_{p, f}}{\sigma_f}\bigg).
\end{equation}

\subsection{Weight calculations}
The weights are calculated by considering a change to the underlying cross sections and then determining whether the event was \textit{more} or \textit{less} probable under that change. Geant4Reweight has the ability to reweight both the inelastic and elastic hadronic cross sections ($\sigma_I$ and $\sigma_E$ respectively).
We consider general variations to the inelastic and elastic cross sections as such:
\begin{equation}
\begin{split}
    \sigma_I &\rightarrow \sigma_I' = f_I\sigma_I
    \\
    \sigma_E &\rightarrow \sigma_E' = f_E\sigma_E.
\end{split}
\end{equation}
A weight is given to an entire track by considering each step, creating a weight for each step, and then multiplying these together:
\begin{equation}
W_{Track} = \prod\limits_i^{f}W_i.
\end{equation}
For each step, one of four things can happen:
\begin{enumerate}
    \item The particle undergoes no interaction
    \item The particle undergoes an inelastic interaction
    \item The particle undergoes an elastic interaction
    \item The particle undergoes some other interaction
\end{enumerate}
Each of these cases have different forms of weights. For the first, the weight is given by the change in the survival probability for the step:
\begin{equation}
    \begin{split}
    W_{S, i} &= e^{-\Delta L_i(\sigma_{I, i}' + \sigma_{E, i}' + \sigma_{X, i})}/e^{-\Delta L_i(\sigma_{I, i} + \sigma_{E, i} + \sigma_{X, i})}
    \\
    &= e^{-\Delta L_i(\sigma_{I, i}' + \sigma_{E, i}')}/e^{-\Delta L_i(\sigma_{I, i} + \sigma_{E, i})}
    \end{split}.
\end{equation}
The other three weights are given by changes to interacting probabilities for that step, as in the last two terms in equation \ref{eqn:p_interact_spec}.
\begin{equation}
    \begin{split}
    W_{I/E, i} &= \bigg(\frac{\sigma_{I/E, i}'}{\sigma_i'}\bigg)
                 \bigg(\frac{\sigma_i}{\sigma_{I/E, i}}\bigg)
                 \bigg(\frac{1 - e^{-\Delta L_i\sigma_i'}}{1 - e^{-\Delta L_i\sigma_i}}\bigg)
                 \\
              &= f_{I/E}\bigg(\frac{\sigma_i}{\sigma_i'}\bigg)
                 \bigg(\frac{1 - e^{-\Delta L_i\sigma_i'}}{1 - e^{-\Delta L_i\sigma_i}}\bigg)
    \end{split}
\end{equation}
\begin{equation}
    \begin{split}
    W_{X, i} &= \bigg(\frac{\sigma_{X, i}}{\sigma_i'}\bigg)
             \bigg(\frac{\sigma_i}{\sigma_{X, i}}\bigg)
             \bigg(\frac{1 - e^{-\Delta L_i\sigma_i'}}{1 - e^{-\Delta L_i\sigma_i}}\bigg)
             \\
          &= \bigg(\frac{\sigma_i}{\sigma_i'}\bigg)
              \bigg(\frac{1 - e^{-\Delta L_i\sigma_i'}}{1 - e^{-\Delta L_i\sigma_i}}\bigg)
    \end{split}
\end{equation}
Here, $\sigma_i = \sigma_{I, i} + \sigma_{E, i} + \sigma_{X, i}$ and $\sigma_i' = \sigma_{I, i}' + \sigma_{E, i}' + \sigma_{X, i}$.

\subsection{Exclusive channels}\label{subsec:exc_chan}
In neutrino experiments it is often important to consider not only the total (also called inclusive) inelastic hadronic interactions but also exclusive interactions. The examples listed in section \ref{sec:intro} highlight this. For examples 1 and 2, the instigating interaction is $\pi^+$ absorption. In example 3, it is $\pi^+$ charge exchange. However, in Geant4 the concept of exclusive hadronic interactions is not explicitly defined, and the only accessible information is the outgoing particles and their kinematics after a model is invoked. In this way, the model can be thought of as a ``black box'' that determines the results of the interaction (outgoing particles) which are then added to the simulation stack and whose information users are able to access. Therefore, Geant4Reweight defines exclusive interaction channels by the outgoing particles following an interaction. For pions, the exclusive channel names and their definitions are given in table \ref{tab:pi_exc}.

\begin{table}[htbp]
\centering
\caption{\label{tab:pi_exc} Summary of exclusive channel definitions for charged pions. $X$ represents any number of nucleons in the final state. $n\pi$ represents more than one $\pi$ (any charge) in the final state.}
\smallskip
\begin{tabular}{|l|l|}
\hline
Channel Name & Definition \\
\hline
Quasielastic & $\pi^{\pm} + N \rightarrow \pi^{\pm} + N' + X$\\
Absorption & $\pi^{\pm} + N \rightarrow N' + X$\\
Charge Exchange & $\pi^{\pm} + N \rightarrow \pi^0 + N' + X$\\
Double Charge Exchange & $\pi^{\pm} + N \rightarrow \pi^{\mp} + N' + X$\\
Production & $\pi^{\pm} + N \rightarrow n\pi + N' + X$\\
\hline
\end{tabular}
\end{table}

For neutrino oscillation experiments such as those listed previously, the model of choice is the Bertini Cascade model. In this model, the incident hadron is injected into the nucleus and allowed to take a series of steps within the nuclear medium. At each step, the particle has a chance to interact with a particle within the nucleus and add it to the stack of active particles. These then step through the nuclear matter, possibly interacting along the way and activating more particles. This so-called cascade of particles continues until all active particles stop within or exit the nucleus and are added to the transport simulation outside of the nucleus. The latter are the interaction results that users can access and which Geant4Reweight uses to define the exculsive interaction channel. Geant4Reweight includes an application which produces effective exclusive cross sections from the Bertini Cascade that are used to calculate weights for exclusive interactions. It does this by invoking the Bertini Cascade model many times over a given momentum range to determine the fraction of times it ends in a specific final state. This fraction is then multiplied by the inelastic cross section to create exclusive cross sections. This is demonstrated in figure \ref{fig:fracs}, where the exclusive channel fractions are shown on the left and the exclusive cross sections are shown on the right along with the total inelastic cross section. These effective exclusive cross sections can be varied and reweighted in a similar way to the full inelastic cross section.

\begin{figure}[htbp]
    \centering
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/fractions_1_14_21.pdf}
        \caption{Exclusive final state fractions.}
    \end{subfigure}
    ~
    \begin{subfigure}[t]{0.45\textwidth}
        \centering
        \includegraphics[width=\textwidth]{Figures/lar_exclusive_xsec_1_14_21.pdf}
        \caption{Exclusive \& Total Inelastic cross sections}
    \end{subfigure}
    \caption{Bertini Cascade model predictions and Geant4 cross section predictions for inelastic $\pi^+$-Ar interactions. The exclusive cross sections on the right are produced by multiplying the fractions on the left by the total inelastic cross section on the right.}
    \label{fig:fracs}
\end{figure}

The reweighting scheme is easily extended to include the exclusive cross sections. We treat the inelastic cross section as the sum over multiple exclusive channels:
\begin{equation}
    \sigma_I = \sum\limits_c^{N_c}\sigma_c.
\end{equation}
We allow users to specify a variation factor to each of the exclusive channels as well as to the total inelastic cross section similar to before:
\begin{equation}
    \sigma_I \rightarrow \sigma_I' = f_I\bigg(\sum\limits_c^{N_c}f_c\sigma_c\bigg).
\end{equation}
$W_S$, $W_E$, and $W_X$ remain largely the same with $\sigma'$ and $\sigma_{I}'$ changed accordingly. $W_I$, however, is modified to account for the separate variations to the exclusive channels. Here, we also have to account for the exclusive channel ($j$) of the interaction. For clarity, the subscript $i$ denoting that this occurs for a given step has been omitted.
\begin{equation}
    \begin{split}
    W_{I,j} &= \bigg(\frac{f_j \sigma_j}{\sum\limits_c^{N_c} f_c\sigma_c}\bigg)
                 \bigg(\frac{\sigma_I}{\sigma_j}\bigg)
                 \bigg(\frac{f_I\sum\limits_c^{N_c} f_c \sigma_c}{\sigma'}\bigg)
                 \bigg(\frac{\sigma}{\sigma_I}\bigg)
                 \bigg(\frac{1 - e^{-\Delta L\sigma'}}{1 - e^{-\Delta L\sigma}}\bigg)
    \\
          &= f_j f_I \bigg(\frac{\sigma}{\sigma'}\bigg)
             \bigg(\frac{1 - e^{-\Delta L\sigma'}}{1 - e^{-\Delta L\sigma}}\bigg)
    \end{split}
\end{equation}