A Modified Slicing Method with Multi-Dimensional Unfolding to Measure Hadron-Argon Cross Sections

1. Introduction

Liquid argon (LAr) technology is widely used by many previous and current neutrino experiments, such as MicroBooNE [1], ArgoNeuT [2], ICARUS [3], and it is also planned to be employed by future experiments such as SBND [4], as well as one of the next-generation large-scale neutrino experiments DUNE [5]. Neutrinos are mostly detected by their interactions with argon nuclei, in which many types of hadrons are involved, including both the knockout of nucleons, and the production of mesons. These neutrino-induced hadrons can also interact with nucleons before they escape the nucleus. This may change the kinematics and types of final-state particles that are detected, which complicates the reconstruction of neutrino interactions. These are known as final-state interaction (FSI) effects. In addition, the propagation and further interactions of these final-state hadrons also need to be simulated properly. Therefore, the knowledge of hadron-argon cross sections is required, which is useful to inform FSI and improve simulations as well as their associated uncertainties.

However, there is little experimental data available on argon, and the predictions are mostly derived by interpolating cross-section results from lighter and heavier nuclei [6], such as carbon, sulfur, and iron, which have more data available [7,8,9,10]. In those experiments, the common setup is to implement a beam of a certain type of hadron of interest, and shoot the beam towards a thin target of the material of interest. The survival rate of the beam hadron after the thin target can be measured and be used to calculate the cross section. The increasing popularity of LAr-based detectors has motivated efforts towards making cross-section measurements on LAr. The LArIAT collaboration proposed the thin-slice method to measure hadron-argon cross section using a LAr time projection chamber (LArTPC) [11], which itself can no longer be considered as a thin target of LAr for hadrons. The precise track reconstruction of LArTPC enables researchers to hypothetically divide the detector into several thin slices, and each slice can be considered as an independent thin-target measurement.

The original method treats the measured cross section in each bin independently, and performs an effective correction in each bin to account for inefficiency and bin migration caused by detector resolution. We keep the essential idea of the thin-slice method, and further develop the method with more rigorous statistical procedures, including using multi-dimensional unfolding to consider the full correlations of different measurements. In this paper, Section 2 shows the derivation of the cross-section formula, and Section 3 describes the slicing method in more details. Section 4 describes the measurement procedures on a simplified simulation sample, where all results are derived using an IPython notebook, referred to as hadron-Ar_XS [12]. Some further discussions on the results as well as summary are given in Section 5.

2. Cross-section formula

The total cross section ${\sigma }_{\mathrm{total}}$ as a function of the incident particle’s kinetic energy E1 is defined according to

$$\frac{d\mathsf{\Phi }}{dx}=-n{\sigma }_{\mathrm{total}}\mathsf{\Phi },$$

(1)

where $\mathsf{\Phi }$ denotes the particle beam flux, and $d\mathsf{\Phi }$ is the infinitesimal reduction of flux. n is the number density of the target material. By moving $dx$, the infinitesimal path length of the particle inside the material, to the right-hand side of Equation 1, and then integrating both sides, we get

$$\mathsf{\Phi }={\mathsf{\Phi }}_0{e}^{-n{\sigma }_{\mathrm{total}}\delta x},$$

(2)

where $\delta x$ is the path length integral. This assumes the cross section $\sigma$ remains constant within the variation of E during its passage of $\delta x$.2 For a certain area and a certain period of time, the number of surviving particles detected is proportional to the outgoing particle flux, and thus we have

$$e^{-n{\sigma }_{\mathrm{total}}\delta x}=\frac{\mathsf{\Phi }}{{\mathsf{\Phi }}_0}=\frac{N_{\mathrm{surviving}}}{N_{\mathrm{incident}}}.$$

(3)

We can also define the number of interacting particles as $N_{\mathrm{interacting}}=N_{\mathrm{incident}}-N_{\mathrm{surviving}}$. Therefore, after measuring the number of incident particles and the number of surviving particles, the total cross section can be calculated as

$${\sigma }_{\mathrm{total}}=-\frac{1}{n\delta x}ln\left(\frac{N_{\mathrm{surviving}}}{N_{\mathrm{incident}}}\right)=-\frac{1}{n\delta x}ln\left(\frac{N_{\mathrm{incident}}-N_{\mathrm{interacting}}}{N_{\mathrm{incident}}}\right).$$

(4)

When it comes to the exclusive cross section3, we denote a as the signal interaction, and b as all the other interactions, and thus we have $d\mathsf{\Phi }=d{\mathsf{\Phi }}_a+d{\mathsf{\Phi }}_b$, where $d{\mathsf{\Phi }}_a$ is the reduction of flux due to the signal interaction. Also we have ${\sigma }_{\mathrm{total}}={\sigma }_a+{\sigma }_b$. Separating $N_{\mathrm{interacting}}$ based on the type of interactions into $N_{\mathrm{interacting}}^a+N_{\mathrm{interacting}}^b$ in Equation 4, we get

$${\sigma }_a+{\sigma }_b=-\frac{1}{n\delta x}ln\left(1-\frac{N_{\mathrm{interacting}}}{N_{\mathrm{incident}}}\right)=-\frac{1}{n\delta x}ln\left(1-\frac{N_{\mathrm{interacting}}^a+N_{\mathrm{interacting}}^b}{N_{\mathrm{incident}}}\right).$$

(5)

${\sigma }_a$ and ${\sigma }_b$ in Equation 5 are not separable given the logarithm function at the right-hand side. Only when $N_{\mathrm{interacting}}\ll N_{\mathrm{incident}}$, which implies that $\delta x$ is very small, and the thin-target approximation holds, then we can use the approximation ${lim}_{x\to 0}ln(1-x)=-x$ and get

$${\sigma }_a+{\sigma }_b=\frac{1}{n\delta x}\frac{N_{\mathrm{interacting}}^a+N_{\mathrm{interacting}}^b}{N_{\mathrm{incident}}}.$$

(6)

Therefore, we have

$${\sigma }_a=\frac{1}{n\delta x}\frac{N_{\mathrm{interacting}}^a}{N_{\mathrm{incident}}}(\mathrm{and}{\sigma }_b=\frac{1}{n\delta x}\frac{N_{\mathrm{interacting}}^b}{N_{\mathrm{incident}}}),$$

(7)

which is in fact a direct implication from the definition of exclusive cross-section

$$\frac{d{\mathsf{\Phi }}_a}{dx}=-n{\sigma }_a\mathsf{\Phi }.$$

(8)

However, in the slicing method described in Section 3, $\delta x$ in each slice we used to calculate $\sigma$ is not necessarily small. Therefore, we seek to get an unbiased cross-section formula without the thin-target approximation. From Equation 1 and 8, we have

$$\frac{{\sigma }_{\mathrm{total}}}{{\sigma }_a}=\frac{d\mathsf{\Phi }}{d{\mathsf{\Phi }}_a}.$$

(9)

For a finite $\delta x$, we can estimate this relationship as

$$\frac{{\overline{\sigma }}_{\mathrm{total}}}{{\overline{\sigma }}_a}=\frac{\int d\mathsf{\Phi }}{\int d{\mathsf{\Phi }}_a}=\frac{\Delta \mathsf{\Phi }}{\Delta {\mathsf{\Phi }}_a}=\frac{N_{\mathrm{interacting}}}{N_{\mathrm{interacting}}^a},$$

(10)

where $\overline{\sigma }$ is the effective mean value for the cross section within the variation of E during the passage of $\delta x$.4 Therefore, combining with Equation 4, we have the expression for any channel a as

$${\overline{\sigma }}_a=\frac{N_{\mathrm{interacting}}^a}{n{N}_{\mathrm{interacting}}\delta x}ln\left(\frac{N_{\mathrm{incident}}}{N_{\mathrm{incident}}-N_{\mathrm{interacting}}}\right).$$

(11)

Because we can never measure $\sigma$ in an infinitely small E bin, we will express $\overline{\sigma }$ as $\sigma$ in the following sections. In the thin-target approximation, $N_{\mathrm{interacting}}\ll N_{\mathrm{incident}}$, Equation 11 can be approximated to Equation 7.

3. Slicing method

A LArTPC cannot be seen as a thin target in terms of hadrons, whose mean free path in LAr is on the order of 10 to 100 cm. However, thanks to its high-resolution track reconstruction ability, the LArIAT collaboration proposed the thin-slice method [11], where they hypothetically divide the detector into several slices along the hadron beam direction. Each slice is viewed as a thin target, with a width of several millimeters based on the spacing of the sensing wires. When detecting tracks in the TPC, each slice serves as an independent thin-target experiment. By detecting where the track ends, we know where the interaction happens, and thus fill in the corresponding energy bins of $N_{\mathrm{interacting}}$ and $N_{\mathrm{incident}}$, which are used to calculate the cross section. The final results are rebinned to wider energy bins such as 50 MeV in order to gain statistics.

Based on the thin-slice method, Ref [13] based on a study on the ProtoDUNE-SP experiment [14] first proposed the idea of energy-slicing, where each energy bin is directly considered as a slice, since the cross section is measured as a function of the kinetic energy of the incident particle.5 Figure 1 shows an illustration of a LArTPC. A beam hadron is incident from the left side of the detector, and leaves a track inside the detector. The beam hadron track ends at the end vertex, where either an interaction occurs or the hadron comes to rest, potentially producing some daughter particles, which can be used to determine the type of the interaction. The kinetic energy of the beam hadron when it enters the detector is denoted as $E_{\mathrm{initial}}$, which is known from beam and it approximately follows a Gaussian distribution given the momentum spread. The kinetic energy of the beam hadron at the end vertex is denoted as $E_{\mathrm{end}}$. Given these two energies, the track can be divided into several slices based on the pre-defined energy bins. The bin edges are indicated by dark red bars in Figure 1, where the last bar is dashed because the beam hadron does not reach that energy. As shown in Figure 1, the first complete slice is referred to as the initial slice; the slice which has the end vertex is referred to as the end slice. If the interaction occurring at the end vertex is a signal interaction, then the end slice is also referred to as the interaction slice.

The piece of track prior to the initial slice is referred to as an incomplete slice, which will not be used. On the contrary, $E_{\mathrm{end}}$ is inside the end slice.

For convenience, we define the slice index $\mathrm{ID}$ from 1 to the number of energy bins N, starting with the highest energy bin. Therefore, for each beam hadron track, there is an initial slice index ${\mathrm{ID}}_{\mathrm{ini}}$, an end slice index ${\mathrm{ID}}_{\mathrm{end}}$, as well as an interaction slice index ${\mathrm{ID}}_{\mathrm{int}}$, which is assigned as null if the interaction occurring at the end vertex is not the signal interaction. In addition, if the end vertex is inside the incomplete slice, then the whole track is not usable, and thus the indices for all three slices will be assigned as null. For a sample of events with a beam hadron track in the detector, the distribution of ${\mathrm{ID}}_{\mathrm{ini}}$ forms the initial histogram $N_{\mathrm{ini}}\left(\mathrm{ID}\right)$, and similarly we have the end histogram $N_{\mathrm{end}}\left(\mathrm{ID}\right)$ and the interaction histogram $N_{\mathrm{int}}\left(\mathrm{ID}\right)$. We also define the incident histogram $N_{\mathrm{inc}}\left(\mathrm{ID}\right)$, which will later appear in the cross section formula 14. Each bin of $N_{\mathrm{inc}}\left(\mathrm{ID}\right)$ counts the number of tracks which reach the energy corresponding to the slice index $\mathrm{ID}$, and thus we say the tracks are incident to that slice. Note, for $N_{\mathrm{inc}}\left(\mathrm{ID}\right)$, one event is likely to contribute to multiple bins, since a track can be incident to a sequence of slices until it interacts. In the thin-slice method, for each track, we fill ${\mathrm{ID}}_{\mathrm{ini}}$ into $N_{\mathrm{ini}}\left(\mathrm{ID}\right)$ and ${\mathrm{ID}}_{\mathrm{end}}$ into $N_{\mathrm{end}}\left(\mathrm{ID}\right)$, while $N_{\mathrm{inc}}\left(\mathrm{ID}\right)$ is filled from ${\mathrm{ID}}_{\mathrm{ini}}$ to ${\mathrm{ID}}_{\mathrm{end}}$. $N_{\mathrm{inc}}\left(\mathrm{ID}\right)$ can also be calculated by the derived $N_{\mathrm{ini}}\left(\mathrm{ID}\right)$ and $N_{\mathrm{end}}\left(\mathrm{ID}\right)$ as 6 [13]

$$\begin{array}{c}\hfill N_{\mathrm{inc}}\left(\mathrm{ID}\right)=\sum _{j=\mathrm{ID}}^N{N}_{\mathrm{end}}\left(j\right)-\sum _{j=\mathrm{ID}+1}^N{N}_{\mathrm{ini}}\left(j\right),\\ \hfill \mathrm{or}{N}_{\mathrm{inc}}\left(\mathrm{ID}\right)=\sum _{j=1}^{\mathrm{ID}}{N}_{\mathrm{ini}}\left(j\right)-\sum _{j=1}^{\mathrm{ID}-1}{N}_{\mathrm{end}}\left(j\right).\end{array}$$

(12)

The two expressions are equivalent given

$$\sum _{\mathrm{ID}=1}^N{N}_{\mathrm{ini}}\left(\mathrm{ID}\right)=\sum _{\mathrm{ID}=1}^N{N}_{\mathrm{end}}\left(\mathrm{ID}\right),$$

(13)

which equals the total number of beam hadron events. Given the relationship between the slice index $\mathrm{ID}$ and the energy E by definition, all of these histograms can also be given as energy histograms.

Comparing to Equation 11, replacing $N_{\mathrm{incident}}$ with $N_{\mathrm{inc}}\left(E\right)$, $N_{\mathrm{interacting}}$ with $N_{\mathrm{end}}\left(E\right)$, $N_{\mathrm{interacting}}^a$ with $N_{\mathrm{int}}\left(E\right)$, and also, $\frac{1}{\delta x}$ with $\frac{1}{\delta E}\frac{dE}{dx}$, we derive the cross section for the signal interaction in each energy bin is given by

$$\sigma \left(E\right)=\frac{N_{\mathrm{int}}\left(E\right)}{n{N}_{\mathrm{end}}\left(E\right)\delta E}\frac{dE}{dx}\left(E\right)ln\left(\frac{N_{\mathrm{inc}}\left(E\right)}{N_{\mathrm{inc}}\left(E\right)-N_{\mathrm{end}}\left(E\right)}\right),$$

(14)

where $\delta E$ is the energy bin width, and $\frac{dE}{dx}\left(E\right)$ is the stopping power of the hadron in LAr. Therefore, for each beam hadron event, three properties are needed, which are $E_{\mathrm{initial}}$, $E_{\mathrm{end}}$, and whether or not it is the signal interaction, in order to derive the slice index for $N_{\mathrm{ini}}$, $N_{\mathrm{end}}$, and $N_{\mathrm{int}}$. This allows us to treat the three indices as a combined 3D variable, and thus enabling the multi-dimensional unfolding discussed in Section 4.3 and Section 4.4.

4. Procedures and results

4.1. Simulations

All results presented in this paper are obtained from data simulated in hadron-Ar_XS [12]. Although this paper focuses on the method to calculate the cross section, it is worth describing how the simulation is done. Smooth and positive curves are created to represent the hadron-argon cross sections $\sigma$ as functions of the hadron’s kinetic energy E, as shown in Figure 2. The signal cross section ${\sigma }_{\mathrm{sig}}\left(E\right)$ is the cross section for an exclusive channel that we want to measure, and other cross sections ${\sigma }_{\mathrm{oth}}\left(E\right)$ accounts for all the others. The total cross section is given as ${\sigma }_{\mathrm{tot}}\left(E\right)={\sigma }_{\mathrm{sig}}\left(E\right)+{\sigma }_{\mathrm{oth}}\left(E\right)$.

As described in Section 3, for each beam hadron event, we need its initial kinetic energy $E_{\mathrm{ini}}$, the kinetic energy at the end vertex $E_{\mathrm{end}}$, as well as the type of interaction occurring at the end vertex, in order to use the slicing method to measure the cross section. Therefore, in our simplified simulation, we only aim to generate these three properties for each event. For $E_{\mathrm{ini}}$, we generate a random value following a Gaussian distribution for each event, in order to mimic the momentum spread. For the latter two, we simulate the hadron’s passage inside LAr with a customized step size $\Delta x$.7 This means, in each step, we generate a random indicator and decide whether the signal interaction or other interactions happen. If not, it proceeds to the next step, until the hadron reacts or its kinetic energy reaches zero. Thus, it also includes the simulation of particle energy loss inside LAr. We use the Bethe-Bloch formula [16] to model the mean $dE/dx$ curve as a function of hadron kinetic energy, and in each step, we generate a random value following a simplified version of the Landau-Vavilov distribution [16] as the $dE/dx$ value to be used in the step. As a result, the energy loss in each step would be $\Delta E=dE/dx·\Delta x$ if no interaction occurs. The mean $dE/dx$ curve used in the simulation and an example $dE/dx$ distribution derived at $E=400$ MeV where the mean $dE/dx$ approximates 2.10 MeV/cm are shown in Figure 3.

A simulation sample of 100000 events was generated. In this simplified simulation, each event has only three properties relevant to cross-section calculations, which are $E_{\mathrm{ini}}$, $E_{\mathrm{end}}$, and a flag indicating the fate of the hadron. The distributions of these three properties for the simulation sample are shown in Figure 4.8

4.2. Extracting the true cross section

From the three properties associated with each event, we use an even binning with $\delta E=50$ MeV9 and obtain the relevant histograms as described in Section 3: $N_{\mathrm{initial}}$ as the distribution of $E_{\mathrm{ini}}$, $N_{\mathrm{end}}$ as the distribution of $E_{\mathrm{end}}$, $N_{{\mathrm{interaction}}_{\mathrm{ex}}}$ as the distribution of $E_{\mathrm{end}}$ but only for events having the signal interaction. After that, we calculate $N_{\mathrm{incident}}$ using Equation 12. The obtained histograms are shown in Figure 5.

Inserting these histograms into Equation 14, we derive the signal cross section ${\sigma }_{\mathrm{sig}}\left(E\right)$10, as shown in Figure 6. This cross section is calculated originally using the true values of the three properties of each event, and its consistency with the simulation curve suggests the feasibility of the slicing method. For a quantitative comparison of the extracted cross section from the simulation sample against the input curve, we can calculate ${\chi }^2$ given as

$${\chi }^2=\left[\sigma -{\sigma }_{\mathrm{curve}}\right]·V_{\sigma }^{-1}·{\left[\sigma -{\sigma }_{\mathrm{curve}}\right]}^T,$$

(15)

where ${\sigma }_{\mathrm{curve}}$ is a vector of the input cross section evaluated at the middle point in each energy bin, and $V_{\sigma }$ denotes the covariance matrix for the calculated true cross sections $\sigma$. The derivation of $V_{\sigma }$ is described later in Section 4.3. If the reduced ${\chi }^2$, defined as the ${\chi }^2$ value divided by the number of bins $N_{\mathrm{bins}}$, approximates 1, it suggests good consistency of the sample with the curve. In Section 5, toy studies are performed to further validate the slicing method.

4.3. Deriving the statistical uncertainty

In the last subSection 4.2, we described how to extract the true cross section, but not yet how to calculate its statistical uncertainty. In order to do so, we consider the problem in a higher dimensional space. This is because the three histograms $N_{\mathrm{initial}}$, $N_{\mathrm{end}}$, and $N_{{\mathrm{interaction}}_{\mathrm{ex}}}$ directly derived from each event, are not independent of each other. For example, for different values of $E_{\mathrm{ini}}$, the distributions of $E_{\mathrm{end}}$ are supposed to be different. In order to fully consider the correlations among them, we define a combined variable as

$${\mathrm{ID}}_{\mathrm{com}}={\mathrm{ID}}_{\mathrm{ini}}+(N_{\mathrm{bin}}+1)·{\mathrm{ID}}_{\mathrm{end}}+{(N_{\mathrm{bin}}+1)}^2·{\mathrm{ID}}_{{\mathrm{int}}_{\mathrm{ex}}}.$$

(16)

Here we assign $\mathrm{ID}=0$ for events with null value, as described in Section 3, and thus there is one more bin in addition to the number of energy bins $N_{\mathrm{bin}}$ for each of ${\mathrm{ID}}_{\mathrm{ini}}$, ${\mathrm{ID}}_{\mathrm{end}}$, and ${\mathrm{ID}}_{{\mathrm{int}}_{\mathrm{ex}}}$.11 This definition can be understood as flattening a 3D array into a 1D array. The distribution of this combined variable is shown in Figure 7. Since we have 20 bins for all $N_{\mathrm{initial}}$, $N_{\mathrm{end}}$, and $N_{{\mathrm{interaction}}_{\mathrm{ex}}}$, there are ${(20+1)}^3=9261$ bins for ${\mathrm{ID}}_{\mathrm{com}}$. The entry in each bin is derived by counting (weighted) events independently, and thus a (weighted) Poisson error can be assigned as the statistical uncertainty for each bin content.

After that, we follow standard error propagation by constructing Jacobian matrix defined as $J={(\partial f_i/\partial x_j)}_{m\times n}$, where x denotes the original variable and f denotes the variable to be transformed into, and thus the covariance matrix is propagated by $V_f=J·V_x·J^T$. There are three steps to go from the covariance matrix for the combined variable to the covariance matrix for the cross section, which are also described in the caption of Figure 8:

• Firstly, the combined variable is projected back to the three axes, namely ${\mathrm{ID}}_{\mathrm{ini}}$, ${\mathrm{ID}}_{\mathrm{end}}$, ${\mathrm{ID}}_{{\mathrm{int}}_{\mathrm{ex}}}$12, and the covariance matrix for $({\mathrm{ID}}_{\mathrm{ini}};{\mathrm{ID}}_{\mathrm{end}};{\mathrm{ID}}_{{\mathrm{int}}_{\mathrm{ex}}})$ is derived.
• Secondly, the null value bins with $\mathrm{ID}=0$ for the three variable are ignored, and $N_{\mathrm{inc}}$ calculated by Equation 12 replaces $N_{\mathrm{ini}}$, and the covariance matrix for $(N_{\mathrm{inc}};N_{\mathrm{end}};N_{{\mathrm{int}}_{\mathrm{ex}}})$ is derived.
• Thirdly, the cross section is calculated using the three energy histograms by Equation 14, and its covariance matrix can be derived by calculating the derivatives appearing in the Jacobian matrix.

With these covariance matrices, the error bars in Figure 5 and Figure 6 are obtained. Figure 8 shows the correlation matrices for these four stages, allowing better visualization than covariance matrices. As can be seen in Figure 8 (c), the off-diagonal blocks are not empty, which suggests there are correlations among different histograms. In Figure 8 (d), the correlation matrix for $\sigma \left(E\right)$ is diagonal, which suggests the true cross section in each bin is independent. However, this will not be the case for the measured cross section as shown later in Section 4.5, and thus, considering the problem in a higher dimensional space is necessary for a rigorous uncertainty evaluation. 13

4.4. Measurement effects

In practice, the true values of $E_{\mathrm{ini}}$, $E_{\mathrm{end}}$, and the type of interaction are unknown, but they need to be measured, which may result in some measurement effects, including the detector resolution and inefficiency when measuring the energy, as well as the mis-identification of the type of interaction. $E_{\mathrm{ini}}$ is usually measured with external instruments outside the LArTPC, while $E_{\mathrm{end}}$ is derived from $E_{\mathrm{ini}}$ with the reconstructed track information in the LArTPC. We also rely on random variables to model these effects effectively rather than simulating events from the origins of the effects. Firstly, to simulate the selection process, a score is generated for each event. Based on whether the score is larger or smaller than a threshold, the event is kept or rejected. The score is generated following a Gaussian distribution, whose mean parameter depends on the true $E_{\mathrm{ini}}$ of the event, in order for the efficiency to be non-uniform as a more general case. Secondly, for all events that pass the selection, two random variables following different Gaussian distributions are generated for each event, denoted as $E_1$ and $E_2$. $E_1$ is used to imitate the resolution on measuring $E_{\mathrm{ini}}$, while $E_2$ accounts for the resolution effects involved with energy deposition of the reconstructed track. Therefore, we have $E_{\mathrm{ini}}^{\mathrm{meas}}=E_{\mathrm{ini}}^{\mathrm{true}}+E_1$ and $E_{\mathrm{end}}^{\mathrm{meas}}=E_{\mathrm{end}}^{\mathrm{true}}+E_1+E_2$14. Finally, in order to simulate the mis-identification among three type of interactions, or to say event fates, which are "no interaction", "signal interaction", "other interactions", a $3\times 3$ confusion matrix is defined, where each element indicates the possibility of a true fate being recognized as a measured fate. Random numbers are used in order to decide the measured fate of each event according to the defined confusion matrix. As a result, the simulated measurement effects can be seen in Figure 9.

Therefore, for each event in the simulation sample passing the selection, it has three true properties and three measured properties. With these properties, we are able to determine its index for the combined variable ${\mathrm{ID}}_{\mathrm{com}}$ in both the true space and the measured space, and thus we can use the simulation sample to model the response matrix as well as the efficiency plot for ${\mathrm{ID}}_{\mathrm{com}}$, which would later be the inputs for unfolding [17] as well as efficiency correction. Figure 10 shows the distribution of the measured combined variable ${\mathrm{ID}}_{\mathrm{com}}^{\mathrm{meas}}$, calculated also by Equation 16 but using the measured values.

For real data, the true information is not available, and thus we use unfolding together with the efficiency correction to transform from the measured histogram to the estimated true histogram, which is referred to as the unfolded histogram. However, the number of bins for ${\mathrm{ID}}_{\mathrm{com}}$ in the order of $N_{\mathrm{bins}}^3$ can go to over a thousand. To unfold a histogram with such a large number of bins can be very time-consuming. Fortunately, in our case, despite the large number of bins, the histogram is usually sparse, with many empty bins.15 Therefore, we delete these empty bins in the ${\mathrm{ID}}_{\mathrm{com}}^{\mathrm{true}}$ histogram and the ${\mathrm{ID}}_{\mathrm{com}}^{\mathrm{meas}}$ histogram separately, re-index the remaining bins, and denote the new index as ${\mathrm{ID}}_{\mathrm{rem}}$. A map is created between ${\mathrm{ID}}_{\mathrm{com}}^{\mathrm{true}}$ and ${\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{true}}$, and another map is created between ${\mathrm{ID}}_{\mathrm{com}}^{\mathrm{meas}}$ and ${\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{meas}}$. After that, we build the response matrix as well as the efficiency plot of ${\mathrm{ID}}_{\mathrm{rem}}$, as shown in Figure 11.

Denoting the response matrix as $R_{ij}=P({\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{meas}}=j|{\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{true}}=i)$, and the efficiency as ${ϵ}_i=P\left(\mathrm{events}\mathrm{with}{\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{true}}=\mathrm{i}\mathrm{being}\mathrm{selected}\right)$, we have ${\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{meas};\mathrm{sim}}=R·(ϵ·{\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{true};\mathrm{sim}})$, where ${\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{true}}$ holds for the simulation sample. For a data sample, we first use the map for the measured histogram to derive ${\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{meas}}$ from ${\mathrm{ID}}_{\mathrm{com}}^{\mathrm{meas}}$. Then we rely on an unfolding algorithm of choice to estimate the unfolding matrix, denoted as $R^{-1}$16, and thus the unfolded ${\mathrm{ID}}_{\mathrm{rem}}$ histogram for data is ${\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{unfd};\mathrm{data}}=(R^{-1}·{\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{meas};\mathrm{data}})·{ϵ}^{-1}$17. Finally the map for the true histogram is used to transform ${\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{true}}$ back to ${\mathrm{ID}}_{\mathrm{com}}^{\mathrm{true}}$.18

4.5. Fake date result

A sample of 10000 events are generated using the same procedure as the simulation sample, but its true information is not used in order to mimic the real data. After selection, we have 5000 events in this fake data sample. Figure 12 shows the correlation matrices involved in the error propagation from the measured ${\mathrm{ID}}_{\mathrm{rem}}$ histograms to the final cross section results. Comparing to Figure 8 when extracting the true cross section, there are two extra steps, which are unfolding and mapping back to ${\mathrm{ID}}_{\mathrm{rem}}$.

In Figure 12 (a), the correlation matrix for ${\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{meas}}$ is diagonal because in each bin the events are counted independently. Figure 12 (b) is the correlation matrix for ${\mathrm{ID}}_{\mathrm{rem}}^{\mathrm{unfd}}$, which is derived from (a) by the unfolding algorithm of choice. We use the Python version of RooUnfoldBayes [19] as the unfolding algorithm, which employs the d’Agostini method [20]. The unfolding algorithm is treated as a black box in this paper, where we input the measured ${\mathrm{ID}}_{\mathrm{rem}}$ histograms as well as its covariance matrix, and get the output including the unfolded ${\mathrm{ID}}_{\mathrm{rem}}$ histograms as well as its covariance matrix.19 We proceed with the number of iterations as 4, which is the only parameter in d’Agostini’s method for unfolding, and it is not optimized in this work. After that, ${\mathrm{ID}}_{\mathrm{rem}}$ is converted back to ${\mathrm{ID}}_{\mathrm{com}}$ under the map created using the measured histograms, and then the following steps are the same with the error propagation of the true cross section as described in Section 4.3. Figure 13 shows the cross section measured using the fake data sample, whose correlation matrix is shown as Figure 12(f). The correlation matrix is not diagonal, which means the measured cross section in each energy bin is correlated, and thus for the final cross section result, we need to present both the central values as well as its covariance matrix. We can also calculate ${\chi }^2$ of the measured cross section against the simulation curve by Equation 15, whose value is shown in the legend of Figure 13. The error bars in the lower-energy bins tend to be larger mostly because the statistics is smaller, since most hadrons interact before they reach these low energies. This trend can also be seen in the true cross section in Figure 6 for the same reason.

5. Discussions and summary

In the previous section, we described how to extract the true cross section of the simulation sample, as well as how to measure the cross section of a data sample. ${\chi }^2$ is calculated in both cases, which is used to quantify the consistency against the simulation curve. In order to further test the results, we perform toy studies. 400 simulation samples, referred to as toys, each with a sample of 10000 events, are generated in the same way as what is described in Section 4. The true cross section as well as its covariance matrix is calculated in each toy simulation sample. For each cross-section bin, we calculate the pull value of each toy, which is defined as

$$\mathrm{pull}=\frac{\sigma \left(E\right)-{\sigma }_{\mathrm{curve}}\left(E\right)}{\sqrt{V_{\sigma }(E,E)}},$$

(17)

where $\sqrt{V_{\sigma }(E,E)}$ is the uncertainty for $\sigma \left(E\right)$. In each bin, the pull values are expected to follow a normal distribution. Figure 14 (a) shows the test results, where a Gaussian distribution is fitted on the pull histograms in each cross-section bin as shown as the blue error bars. By visually comparing with the reference lines, we can see they are generally consistent with the expectation that each of them centers at 0 and has a bar length of 1, corresponding to the two parameters in the Gaussian fit. As another test, which takes into account the covariance among cross-section bins, we show in Figure 14 (b) the histogram of ${\chi }^2$ calculated according to Equation 15 for each toy. A ${\chi }^2$ distribution is fitted on the histogram, whose degree of freedom $N_{\mathrm{df}}$, as shown in the legend, is consistent with the expectation, which is 18 as the number of cross-section bins. These tests serve as a validation of the slicing method. They also suggests that given the current statistics of events in each toy, the bias caused by the imperfection of the model20 is insignificant.

Similarly, we generate 400 toy fake data samples, each with a sample of 10000 events before selection, in order to study the performance of the procedures to measure the cross section. The 400 toy simulation samples used above are combined into a total of 4000000 events, in order to model the response matrix as well as the efficiency plot for each toy fake data sample, and thus we can ignore the statistical uncertainty of the simulation sample. The cross section is measured for each toy fake data sample, and we can also derive the pull distributions in each cross-section bin as well as the histogram of ${\chi }^2$, as shown in Figure 15. In subplot (a), we can see the lengths of blue error bars are generally consistent with 1, but some of their central points show a small bias away from 0. This bias is also referred to as unfolding error. The general unfolding results is effectively applying a re-smearing matrix on the truth information [21]. Treating the re-smeared truth as truth introduces an unfolding error, which can be eliminated by publishing the re-smearing matrix. In addition, the unfolding error tends to be smaller when the regularization becomes weaker with a greater number of iterations. As a result, since we do not include biases, the derived ${\chi }^2$ distribution in Figure 15 (b) is larger.

In the fake data toy study, the simulation sample used to model the response matrix and the efficiency plot is consistent with the toy fake data samples, because they are generated in the same way. When it comes to real data, we need to consider the uncertainties caused by the differences between data and simulation, which can be estimated by fluctuating the relevant parameters of the simulation sample. Additional model validation procedure is essential to examine the compatibility between data and simulation to ensure the differences are within the quoted simulation uncertainties.

In summary, a method as well as the corresponding procedures for the hadron-argon cross section measurement in a LArTPC detector is provided in this paper. The method requires the inputs of the initial energy and the energy at the end vertex of the track, as well as whether it is signal interaction occurring at the end vertex. The method shows good statistical performance, with no obvious bias except for that caused by unfolding, and good estimation of statistical uncertainties, as suggested by the toy studies. To apply it to real data, the systematic uncertainties due to the difference between data and simulation should be considered, and the parameters of the unfolding algorithm used should be optimized with further investigations into the trade-off between bias and variance. These features can be added to the IPython notebook hadron-Ar_XS [12], which also has the potential to be extended to more cross-section studies.