A modified slicing method with multi-dimensional unfolding to measure hadron-argon cross sections

Liquid argon technology is widely used by many previous and current neutrino experiments, and it is also promising for future large-scale neutrino experiments. When detecting neutrinos using liquid argon, many hadrons are involved, which can also interact with argon nuclei. In order to gain a better understanding of the detection processes, and to simulate neutrino events, knowledge of hadron-argon cross sections is needed. This paper describes a new procedure which improves upon the previous work with multi-dimensional unfolding to measure hadron-argon cross sections in a liquid argon time projection chamber. Through a simplified version of simulation, we demonstrate the validity of this procedure.


Introduction
Liquid argon (LAr) technology is widely used by many previous and current neutrino experiments, such as MicroBooNE [1], ArgoNeuT [2], and 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 neutrinoinduced 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, knowledge of hadron-argon cross sections is required, which is useful for informing FSI and improve simulations as well as their associated uncertainties.
However, there are few experimental data available on argon, and the predictions are mostly derived by interpolating cross section results from lighter and heavier nuclei [6,7], such as carbon, sulfur, and iron, which have more data available [8-11].In those experiments, the common set-up was to implement a beam of a certain type of hadron of interest and shoot the beam toward a thin target of the material of interest.The survival rate of the hadron beam after the thin target can be measured and used to calculate the cross section.The increasing popularity of LAr-based detectors has motivated efforts toward making cross section measurements of LAr.The LArIAT collaboration proposed the thin slice method to measure hadron-argon cross sections using an LAr time projection chamber (LArTPC) [12], which itself can no longer be considered a thin target of LAr for hadrons.The precise track reconstruction capability of LArTPCs enables researchers to hypothetically divide the detector into several thin slices, and each slice can be considered an individual thin target experiment.The measurements from all these slices fill the corresponding energy bins.
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 the detector's 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 detail.Section 4 describes the measurement procedures for a simplified simulation sample, where all results are derived using an IPython notebook, referred to as hadron-Ar_XS [13].Some further discussions on the results as well as a summary are given in Section 5.

Cross Section Formula
The total cross section σ total as a function of the incident particle's kinetic energy E1 is defined according to where Φ denotes the particle beam flux, dΦ is the infinitesimal reduction of flux, and n is the number density of the target material.By moving dx, the infinitesimal path length of the particle inside the material shown on the right-hand side of Equation ( 1), and then integrating both sides, we obtain where δx is the path length integral.This assumes that the cross section σ remains constant within the variation of E during its passage of δ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 We can also define the number of interacting particles as N interacting = N incident − N surviving .Therefore, after measuring the number of incident particles and the number of surviving particles, the total cross section can be calculated as follows: 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Φ = dΦ a + dΦ b , where dΦ a is the reduction in flux due to the signal interaction.Also, we have σ total = σ a + σ b .By separating N interacting based on the type of interactions into N a interacting + N b interacting in Equation ( 4), we obtain where σ a and σ b in Equation ( 5) are not separable given the logarithm function on the righthand side.Only when N interacting ≪ N incident , which implies that δx is quite small, and the thin target approximation holds can we use the approximation lim x→0 ln(1 − x) = −x and find Therefore, we have which is in fact a direct implication from the definition of the exclusive cross section However, in the slicing method described in Section 3, δx for each slice, which we used to calculate σ, is not necessarily small.Therefore, we wish to obtain an unbiased cross section formula without the thin target approximation.From Equations ( 1) and ( 8), we have For a finite δx, we can estimate this relationship as where σ is the effective mean value for the cross section within the variation of E during the passage of δx. 4 Therefore, combined with Equation (4), we have the expression for any channel a: Because we can never measure σ in an infinitely small E bin, we will express σ as σ in the following sections.In the thin target approximation, where N interacting ≪ N incident , Equation (11) can be approximated to Equation (7).

Slicing Method
An LArTPC cannot be seen as a thin target in terms of hadrons, whose mean free path in LAr is to the order of 10-100 cm.However, thanks to its high-resolution track reconstruction ability, the LArIAT collaboration proposed the thin slice method [12], where the detector is hypothetically divided 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 individual 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 interacting and N 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 obtain the results.
Based on the thin slice method, Stocker et al. first proposed the idea of energy slicing [14] based on a study of the ProtoDUNE-SP experiment [15].In the energy slicing method, each energy bin is directly considered a slice, which is natural since the cross section is measured as a function of the kinetic energy of the incident particle. 5Figure 1 shows an illustration of an 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 interaction or the decay.The kinetic energy of the beam hadron when it enters the detector is denoted by E initial , which is known from the 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 by E end .Given these two energies, the track can be divided into several slices based on the predefined energy bins. 6The bin edges are indicated by dark red bars in Figure 1, where the last bar is dashed because the beam hadron did not reach that energy.As shown in Figure 1, the first complete slice is referred to as the initial slice, and 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.In contrast, E end is inside the end slice.For convenience, we define the slice index 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 ID ini , an end slice index ID end , as well as an interaction slice index ID int , which is designated 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 designated as null.For a sample of events with a beam hadron track in the detector, the distribution of ID ini forms the initial histogram N ini (ID), and similarly, we have the end histogram N end (ID) and the interaction histogram N int (ID).We also define the incident histogram N inc (ID), which will later appear in the cross section formula in Equation (14).Each bin of N inc (ID) counts the number of tracks which reach the energy corresponding to the slice index ID, and thus we say the tracks are incident to that slice.Note that for N inc (ID), 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, ID ini fills N ini (ID), and ID end fills N end (ID), while N inc (ID) should be filled from the value of ID ini to ID end .Equivalently, N inc (ID) can also be calculated by the derived N ini (ID) and N end (ID) as shown below: [14]   7 The two expressions are equivalent given that which equals the total number of beam hadron events.Given the relationship between the slice index ID and the energy E by definition, all of these histograms can also be given as energy histograms.
Compared with Equation ( 11), by replacing N incident with N inc (E), N interacting with N end (E), N a interacting with N int (E), and also 1 δx with 1 δE dE dx , we derive the cross section for the signal interaction in each energy bin, given by where δE is the energy bin width and dE dx (E) is the stopping power of the hadron in LAr.Therefore, for each beam hadron event, three properties are needed, which are E initial , E end , and whether or not there is signal interaction in order to derive the slice indices for N ini , N end , and N int .This allows us to treat the three indices as a combined 3D variable, thus enabling the multi-dimensional unfolding discussed in Sections 4.3 and 4.4.

Procedures and Results
This section describes the detailed procedures for measuring the hadron-argon cross section with the help of the IPython notebook hadron-Ar_XS [13].We first describe how the simulation samples are prepared (Section 4.1) and then use their true information to extract the true cross sections (Section 4.2).The derivation of statistical uncertainty is described in Section 4.3.Section 4.4 talks about how we model the measurement effects and prepare the fake data sample.Finally, the measured cross section results of the fake data sample are presented in Section 4.5.

Simulations
All results presented in this paper were obtained from data simulated in hadron-Ar_XS [13].Although this paper focuses on the method for calculating the cross section, it is worth describing how the simulation was carried out.Smooth and positive curves were created to represent the hadron-argon cross sections σ as functions of the hadron's kinetic energy E, as shown in Figure 2. The signal cross section σ sig (E) was the cross section for an exclusive channel that we wanted to measure, and other cross sections σ oth (E) accounted for all the others.The total cross section is given as σ tot (E) = σ sig (E) + σ oth (E).
As described in Section 3, for each beam hadron event, we needed its initial kinetic energy E ini , the kinetic energy at the end vertex E 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 aimed to generate these three properties for each event.For E ini , we generated a random value following a Gaussian distribution for each event in order to mimic the momentum spread. 8For the latter two, we simulated the hadron's passage inside LAr with a customized step size ∆x. 9This means that in each step, we generated a random indicator and decided whether the signal interaction or other interactions happened.If not, then it proceeded to the next step until the hadron reacted or its kinetic energy reached zero.Thus, it also included the simulation of particle energy loss inside LAr.We used the Bethe-Bloch formula [17] to model the mean dE/dx curve as a function of the hadron kinetic energy.In each step, a random value was generated following a simplified version of the Landau-Vavilov distribution [17] as the dE/dx value to be used in the step.The mean dE/dx value of the Landau-Vavilov distribution employed in each step aligned with the value calculated by the Bethe-Bloch formula at the kinetic energy in that step.As a result, the energy loss in each step would be ∆E = dE/dx • ∆x if no interaction occurred.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 was approximated to be 2.10 MeV/cm, are shown in Figure 3.  Cross section curves based on which the simulation was generated.The total cross section (blue dash-dotted curve) is the sum of the signal cross section (orange solid curve) and the other cross sections (green dashed curve).The curves do not correspond to any real hadrons in LAr, but they were created to have an order of magnitude in the hundreds of millibarns, which is similar to the real case.The signal cross section curve also imitates a ∆ resonant peak at about 200 MeV in the low-energy region.A simulation sample of 100,000 events was generated.In this simplified simulation, each event had only three properties relevant to the cross section calculations, which were E ini , E 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.

Extracting the True Cross Section
From the three properties associated with each event, we used an even binning with δE = 50 MeV 11 and obtained the relevant histograms as described in Section 3: N initial as the distribution of E ini , N end as the distribution of E end , and N interaction ex as the distribution of E end , but only for events having the signal interaction.After that, we calculated N incident using Equation (12).The obtained histograms are shown in Figure 5.
By inserting these histograms into Equation ( 14), we derived the signal cross section σ sig (E) 12 , as shown in Figure 6.This cross section was 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 χ 2 , given as where σ curve is a vector of the input cross section evaluated at the middle point in each energy bin and V σ denotes the covariance matrix for the calculated true cross sections σ.
The derivation of V σ is described later in Section 4.3.The p value shown in the legend in Figure 6 suggests good consistency of the sample with the curve.In Section 5, toy studies are performed to further validate the slicing method.

Deriving the Statistical Uncertainty
In the last Section 4.2, we described how to extract the true cross section but not yet how to calculate its statistical uncertainty.In order to accomplish this, we considered the problem in the multi-dimensional variable space.This is because the three histograms N initial , N end , and N interaction ex directly derived from each event were not independent of each other.For example, for different values of E ini , the distributions of E end are supposed to be different.In order to fully consider the correlations among them, we defined a combined variable: Here, we assigned ID = 0 for events with a null value, which is defined in Section 3, and thus there was one more bin in addition to the number of energy bins N bin for each of ID ini , ID end , and ID int ex . 13This 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 had 20 bins for all of N initial , N end , and N interaction ex , there were (20 + 1) 3 = 9261 bins for ID com .The entry in each bin was derived by counting the (weighted) events independently, and thus a (weighted) Poisson error could be assigned as the statistical uncertainty for each bin content.After that, we followed standard error propagation by constructing a Jacobian matrix, defined J = (∂ f i /∂x j ) m×n , where x denotes the original variable and f denotes the variable it will be transformed into, and thus the covariance matrix is propagated by There are three steps to transition 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 14 :

•
Firstly, the combined variable is projected back to the three axes, namely ID ini , ID end , and ID int ex , 15 and the covariance matrix for (ID ini ; ID end ; ID int ex ) is derived.• Secondly, the null value bins with ID = 0 for the three variables are ignored.N inc , calculated by Equation ( 12), replaces N ini , and the covariance matrix for (N inc ; N end ; N int ex ) is derived.• Thirdly, the cross section is calculated using the three energy histograms with 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 Figures 5 and 6 were obtained.Figure 8 shows the correlation matrices for these four stages, allowing better visualization than the covariance matrices.As can be seen in Figure 8c, the off-diagonal blocks were not empty, which suggests there were correlations among different histograms.In Figure 8d, the correlation matrix for σ(E) is diagonal, which suggests that the true cross section in each bin was independent.However, this would not be the case for the measured cross section, as shown later in Section 4.5.Thus, considering the problem in the multi-dimensional variable space is necessary for a rigorous uncertainty evaluation. 14In Figure 8a, the matrix may look to be empty because it is sparse, and the bins may be too small for readers to visualize its color.It is still kept for consistency with the other subplots.Similarly, this happens to Figure 12c. 15This can be accomplished by calculating

Measurement Effects
In practice, the true values of E ini and E 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 misidentification of the type of interaction.E ini is usually measured with external instruments outside the LArTPC, while E end is derived from E ini with the reconstructed track information in the LArTPC.We also relied 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 was generated for each event.Based on whether the score was larger or smaller than a threshold, the event was kept or rejected.The score was generated following a Gaussian distribution, whose mean parameter depended on the true parameters of the event, in order for the efficiency to be non-uniform as a more general case. 16In all, 40,917 out of 100,000 events in the simulation sample passed selection.Secondly, for all events that passed selection, two random variables following different Gaussian distributions were generated for each event, denoted as E 1 and E 2 .E 1 was used to imitate the resolution for measuring E ini , while E 2 accounted for the resolution effects involved with energy deposition of the reconstructed track.Therefore, we had Finally, in order to simulate the misidentification among the three types of interactions, or to say the event fates, which were "no interaction", "signal interaction", and "other interactions", a 3 × 3 confusion matrix was defined, where each element indicated the possibility of a true fate being recognized as a measured fate.Random numbers were used in order to decide the measured fate of each event according to the defined confusion matrix.As a result, the simulated measurement effects 18 can be seen in Figure 9. Therefore, for each event in the simulation sample passing selection, it had three true properties and three measured properties.With these properties, we were able to determine its index for the combined variable ID com in both the true space and the measured space, and thus we could use the simulation sample to model the response matrix as well as the efficiency plot for ID com , which would later be the inputs for unfolding [20] as well as efficiency correction.Figure 10 shows the distribution of the measured combined variable ID meas com , which was also calculated by Equation ( 16) but with the measured values.For the real data, the true information was not available, and thus we used 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 ID com to the order of N 3 bins could reach over a thousand.To unfold a histogram with such a large number of bins can be quite time-consuming.Fortunately, in our case, despite the large number of bins, the histogram was usually sparse, with many empty bins. 19Therefore, we deleted these empty bins in the ID true com histogram and the ID meas com histogram separately, re-indexed the remaining bins, and denoted the new index as ID rem .A map was created between ID true com and ID true rem , and another map was created between ID meas com and ID meas rem .After that, we built the response matrix as well as the efficiency plot of ID rem , as shown in Figure 11.By denoting the response matrix as R ij = P(ID meas rem = j|ID true rem = i) and the efficiency as ϵ i = P(events with ID true rem = i being selected), we had ID meas;sim rem = R • (ϵ • ID true;sim rem ), where ID true rem held for the simulation sample.For a data sample, we first used the map for the measured histogram to derive ID meas rem from ID meas com .Then, we relied on an unfolding algorithm of choice to estimate the unfolding matrix, denoted as R, 20 and thus the unfolded 19 For example, out of 9261 bins in total, for either the ID true com histogram in Figure 7 or the ID meas com histogram in Figure 10, only 311 bins and 350 bins were non-empty, respectively. 20It might be natural to think of R as the direct inverse of R, but it has been proven to be problematic to use.
This has been described in many references about unfolding, such as [20].
ID rem histogram for the data was ID unfd;data rem 21 Finally, the map for the true histogram was used to transform ID true rem back into ID true com . 22

Fake Data Results
A sample of 10,000 events was generated using the same procedure as the simulation sample, but its true information was not used in order to mimic the real data.After selection, we had 4,190 events in this fake data sample.The selection rate was similar to the simulation sample.Figure 12 shows the correlation matrices involved in the error propagation from the measured ID rem histograms to the final cross section results.Compared with Figure 8, when extracting the true cross section, there were two extra steps, which were unfolding and mapping back to ID rem . 21If ϵ in bin i is zero, then the value in the bin ID unfd;data rem = i can be estimated using the simulation sample normalized to the data sample directly, because the zero efficiency is usually due to low statistics and will not change the final result significantly.However, the uncertainty associated with this can be evaluated by fluctuating these bin entries. 22It is possible that a bin for ID com is empty in the simulation sample but not empty in the data sample, especially for bins with low statistics.In this case, we can add these non-empty bins of data to the map as well.It is not necessary for the map to be the same for all data samples.In Figure 12a, the correlation matrix for ID meas rem is diagonal because in each bin, the events were counted independently.Figure 12b is the correlation matrix for ID unfd rem , which was derived from (a) using the unfolding algorithm of choice.We used the Python version of RooUnfoldBayes [22] as the unfolding algorithm, which employs the d'Agostini method [23].The unfolding algorithm was treated as a black box in this paper, where we input the measured ID rem histograms as well as its covariance matrix and obtained the output, including the unfolded ID rem histograms as well as its covariance matrix. 23e proceeded with the number of iterations being four, which was the only parameter in d'Agostini's method for unfolding, and it was not optimized in this work.After that, ID rem was converted back into ID com under the map created using the measured histograms, and then the following steps were the same as those for 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 in Figure 12f.The correlation matrix is not diagonal, which means the measured cross section in each energy bin was correlated, and thus for the final cross section result, we needed to present both the central values as well as their covariance matrix.We could also calculate χ 2 of the measured cross section against the simulation curve with Equation ( 15), whose value is shown in the legend of Figure 13.The error bars in the lower-energy bins tended to be larger, mostly because the statistics were 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.Figure 13.The cross section measured using the unfolded histograms of the fake data sample.The right-tail p value was calculated assuming a χ 2 distribution with the number of degrees of freedom N df being 18, which is the number of cross section bins.

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.Here, χ 2 was calculated in both cases, which was used to quantify the consistency against the simulation curve.In order to further test the results, we performed toy studies.Four hundred simulation samples, referred to as toys, each with a sample of 10,000 events, were generated in the same way as what was described in Section 4. The true cross section as well as its covariance matrix was calculated in each toy simulation sample.For each cross section bin, we calculated the pull value of each toy, which is defined as where V σ (E, E) is the uncertainty for σ(E).In each bin, the pull values were expected to follow a normal distribution.Figure 14a shows the test results, where a Gaussian distribution was fitted onto the pull histograms in each cross section bin, as shown by the blue error bars.By visually comparison with the reference lines, we can see that they were generally consistent with the expectation that each of them centered at zero and had a bar length of one, corresponding to the two parameters in the Gaussian fit.As another test, which took into account the covariance among cross section bins, we show in Figure 14b the histogram of χ 2 , calculated according to Equation (15), for each toy.A χ 2 distribution was fitted onto the histogram, whose degree of freedom N df , as shown in the legend, was consistent with the expectation, with 18 being the number of cross section bins.These tests served as a validation of the slicing method.They also suggest that given the current statistics of the events for each toy, the bias caused by the approximations of the method24 was insignificant.Similarly, we generated 400 toy fake data samples, each with a sample of 10,000 events before selection, in order to study the performance of the procedures to measure the cross section.The 400 toy simulation samples used above were combined into a total of 4,000,000 events in order to model the response matrix as well as the efficiency plot for each toy fake data sample, and thus we could ignore the statistical uncertainty of the simulation sample.The cross section was measured for each toy fake data sample, and we could also derive the pull distributions in each cross section bin as well as the histogram of χ 2 , as shown in Figure 15.In subplot (a), we can see that the lengths of the blue error bars were generally consistent with one, but some of their central points showed a small bias away from zero.This bias is considered the unfolding error.The general unfolding result effectively applies a re-smearing matrix on the true information [24].Treating the re-smeared truth as the truth introduces an unfolding error, and thus publishing the re-smearing matrix is suggested in order for others to consider this error when comparing the results.In addition, the unfolding error tends to be smaller when the regularization becomes weaker with a greater number of iterations.Since we did not include this bias, the derived χ 2 was supposed to be larger, whose distribution is shown in Figure 15b.
In the fake data toy study, the simulation sample used to model the response matrix and the efficiency plot was consistent with the toy fake data samples because they were generated in the same way.When it comes to real data, we need to consider the uncertainties caused by the differences between the data and simulation, which can be estimated by fluctuating the relevant parameters of the simulation sample.Additional model validation procedures are essential to examine the compatibility between the data and simulation and ensure the differences were within the quoted simulation uncertainties.In summary, a method as well as the corresponding procedures for the hadron-argon cross section measurement in an LArTPC detector were 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 showed 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 the 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 [13], which also has the potential to be extended to more cross section studies.

Figure 1 .
Figure 1.An illustration of an LArTPC, where a beam hadron is shot into the detector from the left side.More descriptions of the elements in the illustration are provided in the text.
mb) User-defined cross sections used in simulation Total cross section Signal cross section Other cross sections

Figure 2 .
Figure2.Cross section curves based on which the simulation was generated.The total cross section (blue dash-dotted curve) is the sum of the signal cross section (orange solid curve) and the other cross sections (green dashed curve).The curves do not correspond to any real hadrons in LAr, but they were created to have an order of magnitude in the hundreds of millibarns, which is similar to the real case.The signal cross section curve also imitates a ∆ resonant peak at about 200 MeV in the low-energy region.

1 Figure 3 .
Figure 3. (a) The mean dE/dx curve used in the simulation.The dashed vertical line at E = 400 MeV indicates the case of the example dE/dx distribution.(b) The example dE/dx distribution at E = 400 MeV, where the mean dE/dx is approximately 2.10 MeV/cm, which aligns with the value given in subplot (a). 10 the flag indicating the fate of the hadron (c)

Figure 4 .
Figure 4.For the simulation sample, the distribution of (a) E ini , (b) E end , and (c) the flag indicating the fate of the hadron, where it either has no interaction before it comes to rest, or it has signal interaction or other interactions.

Figure 5 .Figure 6 .
Figure 5. Energy histograms derived from the simulation sample: (a) N initial (E), (b) N end (E), (c) N interaction ex (E), and (d) N incident (E).The first and last energy bins are given as overflows.The derivation of error bars is described later in Section 4.3.

Figure 7 .
Figure 7.The distribution of the combined variable ID com of the simulation sample.

Figure 8 .
Figure 8.(a) Correlation matrix for the combined variable ID com , which is diagonal since the entry in each bin was derived by counting independently.(b) Correlation matrix for (ID ini ; ID end ; ID int ex ), where the first block (bin indices 0-20) corresponds to ID ini , the second block (bin indices 21-41) corresponds to ID end , and the third block (bin indices 42-62) corresponds to ID int ex .(c) Correlation matrix for (N inc ; N end ; N int ex ), where the first block (bin indices 0-19) corresponds to N inc , the second block (bin indices 20-39) corresponds to N end , and the third block (bin indices 40-59) corresponds to N int ex .(d) Correlation matrix for the extracted cross section σ(E), which has 18 bins on each axis without the underflow and overflow bins.

Figure 9 .
Figure 9. (a) The distribution of E meas ini (orange histogram) and E true ini (blue histogram) for the simulation sample.(b) The distribution of E meas end (orange histogram) and E true end (blue histogram) for the simulation sample.(c) The resulting confusion matrix for the event fates of the simulation sample.The horizontal axis indicates the measured fates, and the vertical axis indicates the true fates.The color bar indicates the (weighted) event counts in each bin, which add up to be the event counts passing selection.

Figure 10 .
Figure 10.The distribution of the measured combined variable ID meas com of the simulation sample.

Figure 11 .
Figure 11.(a) The response matrix modeled using the simulation sample.The horizontal axis is ID meas rem , and the vertical axis is the true ID true rem .The color bar indicates the (weighted) event counts for each bin.(b) The efficiency for each ID true rem bin.The uncertainty for efficiency was calculated according to the Clopper-Pearson method [21].

Figure 12 .
Figure 12.Correlation matrices for (a) ID meas rem , (b) ID unfd rem , (c) ID com , (d) (ID ini ; ID end ; ID int ex ), (e) (N inc ; N end ; N int ex ), and (f) the measured cross section σ(E) for the fake data sample.
using the fake data sample Signal cross section used in simulation Cross section measured using unfolded histograms 2 /N df = 24.58/18(p-value = 0.14)

Figure 14 .
Figure 14.Toy studies using the extracted true cross sections of 400 toy simulation samples.(a) The pull value test results.The horizontal axis is the energy slice index ID, where ID = 2 corresponds to an energy bin of [900, 950] MeV and ID = 19 corresponds to an energy bin of [50, 100] MeV.The vertical axis is the pull value.The green lines are the pull values of each toy in each ID bin.The blue point and its error bars indicate µ and σ of the Gaussian fit in each ID bin.The red lines sandwiching the blue point indicate the fit error of µ in each ID bin, which can be visually compared to the dashed orange line for its consistency with 0. The dark blue lines sandwiching the end points of the blue error bars indicate the fit error of σ in each ID bin, which can be visually compared to the dotted red line for its consistency with 1.(b) The histogram of χ 2 against the simulation curve.Fitted χ 2 distributions using both the maximum likelihood (MLH) fit and the least chi-square (LCS) fit are overlaid, with the result given in the legend.

Figure 15 .
Figure15.Toy studies using the measured cross sections of 400 toy fake data samples.The detailed descriptions of subplots (a,b) are the same as described in Figure14.