Comparison of ν μ − Ar multiplicity distributions observed by MicroBooNE to GENIE model predictions

C. Adams,9 R. An,10 J. Anthony,3 J. Asaadi,27 M. Auger,1 S. Balasubramanian,31 B. Baller,8 C. Barnes,16 G. Barr,19 M. Bass,19, 2 F. Bay,28 A. Bhat,24 K. Bhattacharya,20 M. Bishai,2 A. Blake,12 T. Bolton,11 L. Camilleri,7 D. Caratelli,7 R. Castillo Fernandez,8 F. Cavanna,8 G. Cerati,8 H. Chen,2 Y. Chen,1 E. Church,20 D. Cianci,7 E. Cohen,25 G. H. Collin,15 J. M. Conrad,15 M. Convery,23 L. Cooper-Troendle,31 J. I. Crespo-Anadón,7 M. Del Tutto,19 D. Devitt,12 A. Diaz,15 S. Dytman,21 B. Eberly,23 A. Ereditato,1 L. Escudero Sanchez,3 J. Esquivel,24 J. J Evans,14 A. A. Fadeeva,7 B. T. Fleming,31 W. Foreman,4 A. P. Furmanski,14 D. Garcia-Gamez,14 G. T. Garvey,13 V. Genty,7 D. Goeldi,1 S. Gollapinni,26 E. Gramellini,31 H. Greenlee,8 R. Grosso,5 R. Guenette,9 P. Guzowski,14 A. Hackenburg,31 P. Hamilton,24 O. Hen,15 J. Hewes,14 C. Hill,14 J. Ho,4 G. A. Horton-Smith,11 A. Hourlier,15 E.-C. Huang,13 C. James,8 J. Jan de Vries,3 L. Jiang,21 R. A. Johnson,5 J. Joshi,2 H. Jostlein,8 Y.-J. Jwa,7 D. Kaleko,7 G. Karagiorgi,7 W. Ketchum,8 B. Kirby,2 M. Kirby,8 T. Kobilarcik,8 I. Kreslo,1 Y. Li,2 A. Lister,12 B. R. Littlejohn,10 S. Lockwitz,8 D. Lorca,1 W. C. Louis,13 M. Luethi,1 B. Lundberg,8 X. Luo,31 A. Marchionni,8 S. Marcocci,8 C. Mariani,30 J. Marshall,3 D. A. Martinez Caicedo,10 A. Mastbaum,4 V. Meddage,11 T. Mettler,1 T. Miceli,17 G. B. Mills,13 A. Mogan,26 J. Moon,15 M. Mooney,2, 6 C. D. Moore,8 J. Mousseau,16 M. Murphy,30 R. Murrells,14 D. Naples,21 P. Nienaber,22 J. Nowak,12 O. Palamara,8 V. Pandey,30 V. Paolone,21 A. Papadopoulou,15 V. Papavassiliou,17 S. F. Pate,17 Z. Pavlovic,8 E. Piasetzky,25 D. Porzio,14 G. Pulliam,24 X. Qian,2 J. L. Raaf,8 A. Rafique,11 L. Rochester,23 M. Ross-Lonergan,7 C. Rudolf von Rohr,1 B. Russell,31 D. W. Schmitz,4 A. Schukraft,8 W. Seligman,7 M. H. Shaevitz,7 J. Sinclair,1 A. Smith,3 E. L. Snider,8 M. Soderberg,24 S. Söldner-Rembold,14 S. R. Soleti,19, 9 P. Spentzouris,8 J. Spitz,16 J. St. John,5, 8 T. Strauss,8 K. Sutton,7 S. Sword-Fehlberg,17 A. M. Szelc,14 N. Tagg,18 W. Tang,26 K. Terao,23 M. Thomson,3 M. Toups,8 Y.-T. Tsai,23 S. Tufanli,31 T. Usher,23 W. Van De Pontseele,19, 9 R. G. Van de Water,13 B. Viren,2 M. Weber,1 H. Wei,2 D. A. Wickremasinghe,21 K. Wierman,20 Z. Williams,27 S. Wolbers,8 T. Wongjirad,29 K. Woodruff,17 T. Yang,8 G. Yarbrough,26 L. E. Yates,15 G. P. Zeller,8 J. Zennamo,4 and C. Zhang2

MicroBooNE detector. These data were collected in 2016 with the Fermilab Booster Neutrino Beam, which has an average neutrino energy of 800 MeV, using an exposure corresponding to 5.0 × 10 19 protons-on-target. The analysis employs fully automatic event selection and charged particle track reconstruction and uses a data-driven technique to separate neutrino interactions from cosmic ray background events. We find that GENIE models consistently describe the shapes of a large number of kinematic distributions for fixed observed multiplicity.

I. INTRODUCTION
A growing number of neutrino physics experiments use liquid argon as a neutrino interaction target nucleus in a time projection chamber [1]. Experiments that use or will use liquid argon time projection chamber (LArTPC) technology include those in the Short-Baseline Neutrino (SBN) program [2] at Fermilab, centered on searches for non-standard neutrino oscillations, and the long-baseline DUNE experiment [3]. The SBN program consists of the MicroBooNE experiment [4], an upgraded ICARUS experiment [5], and the new SBND experiment [6]. The DUNE experiment seeks to establish the mass ordering of the three standard model neutrinos and the charge parity violation parameter phase δ CP in the PMNS neutrino mixing matrix [7].
All LArTPC neutrino oscillation-related measurements require a precise understanding of neutrino scattering physics and the measured response of the LArTPC detector to final state particles. These depend on: (a) the neutrino flux seen by the experiment, (b) the neutrino scattering cross sections, (c) the interaction physics of scattering final state particles with argon, (d) transport and instrumentation effects of charge and light in the LArTPCs, and (e) software reconstruction algorithms. In practice item (a) is determined by a combination of hadron production cross section measurements and precise descriptions of neutrino beamline components. Item (b) is most commonly provided by the GENIE [8] neutrino event generation model for neutrino-argon scattering. Items (c)-(e) are incorporated into a detailed suite of GEANT4-based [9] simulation and event reconstruction products called LArSoft [10].
GENIE has been built up from models of the most important physical scattering neutrino-nucleon mechanisms for the SBN and DUNE energy regimes (0.5 − 5 GeV) : quasielastic (QE) scattering ν N → − N , ν N , resonance production (RES) ν N → − R, ν R , and non-resonant multi-hadron production referred to as deep inelastic scattering (DIS): ν N → − X, ν X [11]. The underlying neutrino-nucleon scattering processes receive significant modification from the nuclear environment, including the effects of Fermi motion of target nucleons, many-nucleon effects, and final state interactions (FSI) [12]. While GENIE has received a fair amount of "tuning" (the process of finding a set of GENIE parameters chosen to optimize agreement with a particular data set) from previous electron and neutrino scattering measurements, considerable uncertainties remain in the modeling of both the underlying neutrino-nucleon scattering and the nuclear environment [13].
Relatively few neutrino scattering measurements on argon exist [14][15][16][17][18][19], especially for the recoil hadronic system. Most of these report low-statistics exclusive final states. Nearly all existing neutrino scattering constraints on GENIE models derive measurements on scattering from carbon, which has 30% of argon's atomic mass number and a 22% lower neutron-to-proton ratio. We take a step in improving the empirical understanding of neutrino scattering from argon here by performing a large set of comparisons of observed inclusive properties of charged current scattering, measured at the MicroBooNE experiment in the Fermilab Booster Neutrino Beam (BNB) [20] (MicroBooNE and MiniBooNE share the same beam), to predictions from several variants of GENIE. These comparisons are generated by applying fully automated event reconstruction and signal selection tools to a subset of MicroBooNE's first collected data. While this analysis must focus in large part on reducing cosmic ray backgrounds, sensitivity to GENIE model parameters remains.
In Sec. II, we introduce the formalism that relates interaction models to our measurements. In Sec. III, we describe the MicroBooNE detector and Booster Neutrino Beam. Section IV describes the cosmic ray backgrounds in Micro-BooNE. Section V summarizes the event selection procedure and the main software tools used in the analysis. Section VI presents the procedure for discriminating BNB neutrino interactions from cosmic ray background events. Section VII shows the results of a systematic uncertainty analysis. Section VIII presents the comparison of observed charged particle multiplicity distributions and charged track kinematic distributions for each multiplicity, to predictions from GENIE. Section IX provides a discussion of the result, and Sec. X gives an overall conclusion.

II. INTRODUCTION TO OBSERVED CHARGED PARTICLE MULTIPLICITY AND KINEMATIC DISTRIBUTIONS
Neutrino interactions in the MicroBooNE detector produce charged particles that can be reconstructed as tracks in the liquid argon medium of the MicroBooNE LArTPC. These interactions can be characterized by a number of inclusive properties. The charged particle multiplicity, or number of primary charged particles, n, is a simple observable characterizing final states in high-energy-collision processes, including neutrino interactions. We note that in MicroBooNE the observable charged particle multiplicity corresponds to that of charged particles exiting the target nucleus participating in the neutrino interaction.
The charged particle multiplicity distributions (CPMD) comprise the set of probabilities, P n , associated with producing n charged particles in an event, either in full phase space or in restricted phase space domains. In addition to the observed CPMD, kinematic properties of all charged particle tracks for each multiplicity can be examined. Determination of inclusive event properties such as the CPMD and of individual track kinematic properties at Fermilab BNB neutrino energies naturally fits into the modern strategy [11] of presenting neutrino interaction measurements in the form of directly observable quantities.
Inclusive measurements expand the empirical knowledge of neutrino-argon scattering that will be required by the DUNE experiment and the Fermilab SBN program. As physical observables, the CPMD and other distributions can also be used to test models, or particular tunes of models such as GENIE. These models are typically constructed from a set of exclusive cross section channels, and tests of inclusive distributions can provide independent checks.
We describe here an evaluation of several variants of GE-NIE against observed charged particle distributions, including the observed CPMD in MicroBooNE data collected in 2016 in the Fermilab BNB. For the observed CPMD, we mean the conditional probability, after application of certain detector selection requirements, of observing a neutrino interaction with n charged tracks relative to the probability of observing a neutrino interaction with at least one charged track: where N obs,n is the number of neutrino interaction events with n observed tracks.
Our analysis requires at least one of the charged tracks to be consistent with a muon; hence the O n are effectively observed CPMD for ν µ charged current (ν µ CC) interactions. The ν µ NC, ν e ,ν e , andν µ backgrounds, in total, are expected to be less than 10% of the final sample. The muon candidate is included in the charged particle multiplicity, and all events thus have n ≥ 1. For each multiplicity, we have available the kinematic properties of charged tracks. These can in principle be related to the 4-vector components of each track; however, we choose distributions of directly observable quantities in the detector: visible track length and track angles.
Values for O n depend on cross sections for producing a multiplicity n, σ CC,n , as well as the BNB neutrino flux and detector acceptance and efficiency: . dΠ n dσ CC,n (E ν , Π n ) dΠ n ε n,n (E ν , Π n ) , where E ν is the neutrino energy, Φ ν (E ν ) is the neutrino flux summed over ν µ ,ν µ , ν e , andν e species, dΠ n represents the n -particle final state phase space, ε n,n (E ν , Π n ) is an acceptance and efficiency matrix that gives the probability that an n charged particle final state produced in phase space element dΠ n is observed as an n-particle final state in the detector, and dσ CC,n (E ν , Π n ) /dΠ n are the differential cross sections for producing a multiplicity n . One can likewise express the distribution of any observed kinematic distribution X n corresponding to an observed multiplicity n as whereε n,n (E ν , Π n → X n ) is the probability that an n charged particle final state produced in phase space element dΠ n produces the observed value X n of the kinematic variable in the detector. In practice we obtain the O n and distributions of X n directly from data and compare these to values derived from evaluating Eqs. 2 and 3 using a Monte Carlo simulation that includes GENIE neutrino interaction event generators coupled to detailed GEANT-based models of the Fermilab BNB and the detector. The observed CPMD and inclusive observed kinematic distributions have several desirable attributes. The σ CC,n are all large up to n 4 at these neutrino energies (see Sec. III); therefore only modest event statistics are required. Only minimal kinematic properties of the final state are imposed (the track definition implies an effective minimum kinetic energy), and complexities associated with particle identification and photon reconstruction are avoided. At the same time, the observed quantities reveal much of the power of the LArTPC in identifying and characterizing complex neutrino interactions. The observed CPMD and associated kinematic distribution ratios will have reduced sensitivity to systematic normalization uncertainties associated with flux and efficiency compared to absolute cross section measurements.
A disadvantage of the use of observed CPMD and other kinematic quantities is their lack of portability. One must have access to the full MicroBooNE simulation suite to use the O n to test other models.  . 2) has an active mass of about 85 tons (85 Mg) of liquid argon. It is 10.4 meters long in the beam direction, 2.3 meters tall, and 2.6 meters in the electron drift direction. Electrons require 2.3 ms to drift across the full width of the TPC from cathode (at −70 kV) to anode (at ∼ 0 kV). Events are read out on three anode wire planes with 3 mm spacing between wires. Drifting electrons pass through the first two wire planes, which are oriented at ±60 degrees relative to vertical, producing bipolar induction signals. The third wire plane, the collection plane, has its wires oriented vertically and collects the charge of the drifting electrons in the form of a unipolar signal. The MicroBooNE readout electronics allow for measurement of both the time and charge created by drifting electrons on each wire. The amplified, shaped waveforms from 8256 wires from induction and collection planes are digitized at 2 MHz using 12-bit ADCs. A data acquisition system readout window consisting of 9600 recorded samples (4.8 msec) for all wires is then noise-filtered and deconvolved utilizing offline software algorithms. Reconstruction algorithms are then used on these output waveforms to reconstruct the times and amplitudes of charge depositions (hits) on the wires from particle-induced ionization in the TPC bulk.
While all three anode planes are used for track reconstruction, the collection plane provides the best signal-to-noise performance and charge resolution. The analysis presented here excludes regions of the detector that have non-functional collection plane channels (∼ 10%). It also imposes requirements on the minimum number of collection plane hits−current pulses processed through noise filtering [21], deconvolution, and calibration−associated with the reconstructed tracks. All charged particle track candidates are required to have at least 15 collection plane hits, and the longest muon track candi-date is required to have at least 80 collection plane hits. Furthermore, as described in Sec. V B, we use two discriminants to extract the neutrino interaction and cosmic ray background contributions to our data sample that are based on collection plane hits. A light collection system consisting of 32 8-inch photomultiplier tubes (PMTs) with nanosecond timing resolution enables precise determination of the time of the neutrino interaction, which crucially aids in the reduction of cosmic ray backgrounds.
The BNB employs protons from the Fermilab Booster synchrotron impinging on a beryllium target. The proton beam comes in bunches of protons called "beam spills", has a kinetic energy of 8 GeV, a repetition rate of up to 5 Hz, and is capable of producing 5 × 10 12 protons-per-spill. Secondary pions and kaons decay producing neutrinos with an average energy of 800 MeV. MicroBooNE received 3.6 × 10 20 protons-on-target in its first year of running, from fall 2015 through summer 2016. This analysis uses a fraction of that data corresponding to 5.0 × 10 19 protons-on-target.

IV. COSMIC RAY BACKGROUNDS
The MicroBooNE detector lacks appreciable shielding from cosmic rays (CR) since the detector is at the earth's surface and has little overburden. Most events that are recorded and processed through an online software trigger which requires that the total light recorded with the PMT system exceeds 6.5 photoelectrons (PE) during neutrino beam operations ("on-beam data") contain no neutrino interactions. Triggered events with a neutrino interaction typically have the products of up to 20 cosmic rays coincident with the beam spill in the event readout window (4.8 msec) contributing to a recorded event along with the products of the neutrino collision. A large sample of events recorded under identical conditions as the on-beam data, minus the coincidence requirement with the beam, ("off-beam-data") has been recorded for use in characterizing cosmic ray backgrounds. A straightforward on-beam minus off-beam background subtraction is difficult, as the off-beam data does not reproduce all correlated detector effects associated with on-beam events that contain a neutrino interaction with several overlaid cosmic rays. The situation is particularly complicated with events containing neutrino interactions with N obs,n = 1, which share the same topology with the most common single-muon CR configuration. Monte Carlo simulations of the CR flux using the CORSIKA package [22] provide useful guidance; however, the ability of these simulations to describe the very rare CR topologies that closely match neutrino interactions is not well known.
For these reasons this analysis employs a method to separate neutrino interaction candidates from CR backgrounds that is driven by the data itself. Even though CR tracks should always appear to at least enter the detector, they can satisfy the experimental condition of being fully contained if a segment of the CR track falls outside the data acquisition readout time window, or if a segment of the track fails to be identified due to instrumentation-or algorithm-related inefficiencies. The separation of neutrino interaction candidates from CR backgrounds rests on the observation that a neutrino ν µ CC interaction produces a final state µ − that slows down as it moves away from its production point at the neutrino interaction vertex due to ionization energy loss in the liquid argon. As it slows down, its rate of restricted energy loss [23], dE/dx R , increases, and deviations from a linear trajectory due to multiple Coulomb scattering (MCS) become more pronounced. A CR muon track can produce an apparent neutrino interaction vertex if it comes to rest in the detector or it is not fully reconstructed to the edge of the TPC, but the CR track will exhibit large dE/dx R and MCS effects in the vicinity of this vertex. Furthermore, the vast majority of ν µ CC muons travel in the neutrino beam direction ("upstream" to "downstream"), whereas CR muons move upstream or downstream with equal probability.

A. Data
This analysis uses two data samples: • "On-beam data", taken only during periods when a beam spill from the BNB is actually sent. The on-beam data used in this analysis were recorded from February to April 2016 using data taken in runs in which the BNB and all detector systems functioned well. This sample comprises about 15% of the total neutrino data collected by MicroBooNE in its first running period (October 2015 to summer 2016), • "Off-beam data" taken with the same software trigger settings as the on-beam data, but during periods when no beam was received. The off-beam data were collected from February to October 2016.

B. Simulation
The LArSoft software framework is used for processing data events and Monte Carlo simulation (MC) events in the same way. Three simulation samples are used in this analysis: • Neutrino interactions simulated with a default GENIE model overlaid with CORSIKA CR events ("MC default"), • MC default augmented by the GENIE implementation of the Meson Exchange Current model [24] overlaid with CORSIKA CR events ("MC with MEC"), • MC default augmented by the GENIE implementation of the Transverse Enhancement Model [25] overlaid with CORSIKA CR events ("MC with TEM").
The GENIE models used in this analysis are Relativistic Fermi Gas model [26] for the nuclear momentum, Llewellyn-Smith [27] for QE interactions, Rein-Sehgal model [28] for resonance interactions, and Empirical model by Dytman [29] for MEC interactions. The MC default does not include contributions from the excitation of two particle-two hole (2p2h) [30] final states in neutrino-nucleus scattering, which may be important for low energy neutrinos. The MC with MEC and TEM include the excitation of 2p2h final states in neutrino-nucleus scattering in alternative ways. In TEM, the empirical superscaling function is modeled with an effective spectral function (ESF) [31].
The generator stage (production of a set of final state fourvectors for particles originating from the argon nucleus as a result of the ν µ -Ar interaction in GENIE) employs GE-NIE (version v2.8.6d for the MC default and v2.10.6 for the MC with MEC and TEM) with overlaid simulated CR backgrounds using CORSIKA version v7.4003 with a constant mass composition model (CMC) [32] at 226 m above sealevel elevation. Simulated secondary particle propagation utilizes GEANT version v4.9.6.p04d using a physics particle list QGSP BIC with custom physics list for optical photons, and the detector simulation employs LArSoft version v4.36.00. All GENIE samples were processed with the same GEANT and LArSoft versions for detector simulation and reconstruction. These samples thus allow for relative comparison of different GENIE models to the data.

C. Reconstruction
Event reconstruction makes use of anode plane waveforms from the TPC and light signals from the PMT system. Raw signals from the TPC, recorded on the wires, first pass through a noise filter [21] before hits are extracted from the observed waveforms on each wire. A set of reconstruction algorithms known as Pandora [33,34] combine hits connected in space and time into two-dimensional (2D) clusters for each of the three anode planes. Then a three-dimensional (3D) vertex reconstruction is performed in which the energy deposition around the vertex and the knowledge of the beam direction is used to preferentially select vertex candidates with more deposited energy and with low z coordinates, respectively. 2D clusters in all three planes are then merged into 3D tracks and showers. The preferential direction for tracks and showers is also from the upstream to the downstream end of the detector. In this analysis, for each 3D track candidate, the start position is taken to be the 3D track end closest to the candidate neutrino vertex, if the track satisfies the requirement that the track start position is sufficiently close to the vertex. This analysis makes no use of shower objects, which are mainly associated with electrons and photons.
Light collected on the 32 PMTs in MicroBooNE is used to reconstruct optical hits. To reduce sensitivity to possible fluctuations in the signal baseline, a threshold-based hitreconstruction algorithm requires PMT pulses of a minimum charge for a hit to be reconstructed. A weighted sum of PMT hits (optical flashes) are reconstructed by requiring a time coincidence of ∼1 µs between hits on multiple PMTs. The relative timing and charge collection of optical hits, along with the spatial locations of the PMT, within an optical flash is then used to associate the flash with reconstructed tracks in the TPC, a process known as flash-matching.

D. Event Selection
Event selection starts by requiring an optical flash within the 1.6 µs duration beam spill window and the summed light collected by the PMT to exceed 50 PE. Reconstructed vertices must be contained in the fiducial volume of the detector, defined as 10 cm from the border of the active volume in x, 10 cm from the border of the active volume in z, and 20 cm from the border of the active volume in y (see Fig. 2 for detector coordinates). At each candidate neutrino interaction vertex, a candidate muon track is identified as the longest of all tracks starting within 5 cm of the vertex. The candidate muon track is further required to be fully contained within the detector, where containment requires both ends of the track to lie within the same fiducial volume required for an event vertex, to have at least 75 cm 3D track length, and to have an event vertex located within 80 cm in z of the PMT-reconstructed position of an optical flash. Considerable CR backgrounds remain after these pre-selection procedures, with signal/background ≈ 1/1.
Pre-selected events then pass through a second stage filter that imposes further quality conditions on track candidates. Start and end points of the candidate muon must lie in detector regions with functional collection plane wires. The candidate muon track must start 46 cm below the top surface of the TPC in order to suppress CR backgrounds, must start within 3 cm (reduced from 5 cm) of the selected vertex position, must have at least 80 hits in the collection wire plane, and must not have significant wire gaps in the start and end 20 collection plane-hit segments used in the pulse-height (PH) test (Sec. V E 1) and the multiple Coulomb scattering test (Sec. V E 2).
Events satisfying all of the above criteria constitute the final data sample. Table I lists the event passing rates for the on-beam data, off-beam data, and the MC default samples at different steps of the event selection. The passing rates in on-beam data are consistent with expectations for a mix of CR-only events, as provided by the off-beam data, and events containing a neutrino interaction in addition to cosmic rays, as provided by the MC.
The observed multiplicity of a selected event is defined to be the number of particles starting within 3 cm of the selected vertex that have at least 15 collection plane hits where the Pandora MicroBooNE track reconstruction algorithms perform optimally. There is no containment requirement for tracks other than the candidate muon track. Table II lists the number of selected events in each multiplicity bin with relative event rates for on-beam data, off-beam data, and MC default samples.
The minimum collection plane hit condition corresponds to a minimum range in liquid argon of 4.5 cm, and the requirement thus excludes charged particles below a particle-typedependent kinetic energy threshold from entering our sample that ranges from 31 MeV for a π ± to 69 MeV for a proton. No acceptance exists for particles with kinetic energies below these thresholds, which roughly increase as the secant of the track angle with respect to the neutrino beam direction.
The average efficiency for the Pandora-based track recon-struction used in this analysis is ε ≈ 45% at the 15 collection plane hit threshold. This relatively low value, with implicit kinetic energy thresholds, creates a common occurrence called "feed-down" wherein events produced with n tracks at the argon nucleus exit position are reconstructed with an observed multiplicity n < n. For example, n = 1 is commonly observed because one of the two tracks in a quasi-elastic-like event fails to be reconstructed due to low acceptance or tracking efficiency. The candidate muon containment requirement limits its energy to be 1.2 GeV depending on the muon scattering angle. This results in a sample biased towards relatively higher inelasticity, E H /E ν , with E H being the energy transferred from the neutrino to the hadronic system in the collision. Figure 3 shows the GENIE expectations for the true kinetic energy of muons, protons, and pions produced in BNB neutrino interactions in MicroBooNE. The kinetic energy thresholds associated with the 15 collection plane hit requirement for short tracks and the 75 cm 3D track length requirement for the long track are evident.

E. Event Classification
Selected events are next classified into four categories based on whether they pass or fail the PH test and the MCS test described in the following sections. These are the candidate muon track direction-based tests which are used to separate neutrino signal and CR background contributions in the sample. Table III lists the event selection rates for the on-beam data, off-beam data, and the MC default samples in each category. The final samples are called neutrino-enriched, mixed, or background-enriched sub-samples depending on whether events pass both tests, pass either one of the two tests, or fail both tests, respectively.

Pulse Height Test
A neutrino-induced muon from a CC event will exhibit an increasing rate of energy loss as one moves downstream along its track. A visual diagram for the PH test is shown in Fig. 4. We take into account the expected behavior of the rate of restricted energy loss, dE/dx R , with the following procedure: • Compute the truncated mean of the pulse heights deposited in 20 consecutive collection plane hits, PH U , starting 10 hits away from the upstream end of the muon track that is taken as a proxy for the upstream restricted energy loss. The truncated mean is formed by taking the average of the 20 PH after removing individual PH that do not lie within the range of 20% − 200% of the average [35]: which can be determined iteratively with an initial approximation that PH is the arithmetic average. Use    of the truncated mean PH rather than the average PH minimizes effects of large energy loss fluctuations, • Form a similar quantity from 20 consecutive collection plane hits that end 10 collection plane hits away from the downstream end of the track, PH D , • Form the test p = PH U < PH D . Muons from ν µ CC interactions will pass this test with a probability P (PH).
Muons from CR background can be characterized by the probability that they fail this test, denoted as Q (PH) . tribution for neutrino events only from MC default (signal MC) and off-beam data (cosmic data) samples. The expected signal is considerably enriched relative to the background for PH ratios greater than 1 and we use this value to define the PH test used in the analysis. > 1 pass the PH test.

Multiple Coulomb Scattering Test
A neutrino-induced muon from a CC event will generally exhibit an increasing degree of multiple Coulomb scattering (MCS) about a nominal straight line trajectory as one moves from upstream to downstream along the track. A visual diagram for the MCS test is shown in Fig. 6. We take into account the expected MCS behavior by an independent test with the following procedure: • Compare the hit time predicted at the geometric center of the track, t C , by L M , which uses hits about the geometric center, to the time predicted at the geometric center of the track by the projection of L U from the beginning of the track: • Repeat the process except compare t C from L M to the time predicted at the geometric center of the track by the projection of L D from the end of the track: • Form the test q = ∆t UM < ∆t DM . Since MCS should become, on average, more pronounced along the downstream end of the track as the momentum decreases, this provides a second directional test on the muon track candidate. Muons from ν µ CC interactions will pass this test with a probability P (MCS). Muons from CR background can be characterized by the probability that they fail this test, denoted as Q (MCS) .  We observe that MCS ratio for the signal dominates over the background for values greater than 1 and we use this value to define the MCS test used in this analysis.

A. Data-Driven Signal+Background Model
On-beam data consists of a mixture of neutrino interaction and CR background events. We designate a passing of the PH or MCS test by the symbol ν, and a failure of either of the two tests by the symbol CR, thus creating the categories "νν", "νCR", "CRν", "CRCR", which contain a corresponding number of events N νν , N νCR , N CRν , and N CRCR . The "νν" and "CRCR" categories are expected to have relatively high neutrino or CR purity, respectively; while the "νCR" and "CRν" have mixed purity.
We next build a model for the number of events in each category as follows: The quantitiesN νν ,N CRν ,N νCR , andN CRCR are model parameters corresponding to the observed number of events N νν , N νCR , N CRν , and N CRCR , respectively.N ν andN CR are the estimated number of neutrino and CR events, in the sample, to be determined by a fit described below. The quantities P (PH) and P (MCS) represent the average probabilities that a neutrino interaction muon passes the PH or MCS test condition, while Q (PH) and Q (MCS) denote the mean probabilities that a cosmic ray muon fails one of these tests. The conditional probability P (MCS|PH) denotes the fraction of time that a neutrino interaction muon event that passes the MCS condition after it has passed the PH condition, and the conditional probability Q (MCS|PH) denotes the fraction of time that a cosmic ray event muon fails the MCS test after failing the PH test.
As the MCS and PH conditions result from different physical processes (muon-nucleus and muon-electron scattering, respectively), and the MCS and PH test are formed from different measurements (time and charge, respectively), the PH and MCS tests are nearly independent with P (MCS|PH) ≈ P (MCS) and Q (MCS|PH) ≈ Q (MCS). In the analysis we find evidence for weak, but non-negligible, correlations between the tests, and use the conditional probabilities to take these into account.
We collect data and construct a similar model for off-beam data, which contains no neutrino content, dividing the events into the same categories as above, and fitting the observed number of events in each category, N νν , N νCR , N CRν , and N CRCR to the parameterizations: In this case the νν and CRCR categories are expected to be enriched samples containing muons characteristic of neutrino interactions and cosmic rays, respectively, while the CRν and νCR samples have a mixed composition.N CR is the estimated CR content of the sample (in practice the number of events in the sample).
Our algorithm uses the eight categories of events in onbeam and off-beam data to estimate the neutrino content in each multiplicity bin. To calculate the MC distributions, we replace the on-beam data with the MC samples and perform the fit again. The same off-beam data sample was used in both fits. In the absence of correlations, the quantitiesN ν , N CR ,N CR , P (PH), P (MCS), Q (PH), and Q (MCS) can be directly determined from the data with no model inputs. The addition of the two conditional probabilities P (MCS|PH) and Q (MCS|PH) requires use of a model to determine the correlation between the PH and MCS tests. These correlations are implemented through the parameterizations .
The two new parameters α ν and α CR are obtained from Monte Carlo simulation of neutrino data and from the off-beam data, respectively. If α ν = 1, no correlation would exist between the tests, whereas a large α ν would imply near total correlation, with similar conditions applied to α CR .

B. Fitting Procedure
We construct a likelihood function based on the probability distribution for partitioning events into one of four categories of a multinomial distribution, for both on-beam and off-beam data. The multinomial probability of observing n i events in bin i, with i = 1, 2, 3, 4, with the probability of a single event landing in bin i equal to r i is M (n 1 , n 2 , n 3 , n 4 ; r 1 , r 2 , r 3 , r 4 ) = (n 1 + n 2 + n 3 + n 4 )! n 1 !n 2 !n 3 !n 4 ! r n 1 1 r n 2 2 r n3 3 r n 4 4 .
The n i are the observed number of events in each bin, and the r i are functions of the model parameters.
The likelihood also incorporates the Poisson statistics of observing n 1 + n 2 + n 3 + n 4 in both the on-beam and off-beam data: withN The final likelihood function is The model parametersN ν ,N CR ,N CR , P (PH), P (MCS), Q (PH), and Q (MCS) and their statistical uncertainties are estimated via the maximum likelihood method, implemented by minimizing the negative-log-likelihood using the MIGRAD minimization in the standard MI-NUIT [36] package in ROOT [37].
The fitting procedure can be used to obtain estimates for N ν ,N CR ,N CR , P (PH), P (MCS), Q (PH), and Q (MCS) for each multiplicity. When the probability parameters P (PH), P (MCS), Q (PH), and Q (MCS) are consistent between multiplicities, we use all multiplicities together in their determination for improved statistical precision and vary only the three parametersN ν ,N CR , andN CR for each individual multiplicity.

C. Results with Simulated Events
Maximum likelihood fits were performed on all three GE-NIE simulation samples to extract the values of seven param-etersN ν ,N CR ,N CR , P(PH), Q(PH), P(MCS), and Q(MCS)). Parameters α ν and α CR and their uncertainties were extracted from MC and off-beam data samples and kept fixed for the subsequent fits. As expected, the PH and MCS probabilities show no statistically significant difference between the three GENIE models considered. Table IV lists the values obtained from the fit for the above-mentioned parameters in the default MC and the MicroBooNE data.
The number of neutrino events in the simulated data samples were extracted for each observed multiplicity and compared to the known number from the event generation. Table V and Fig. 8 summarize this comparison. We find that the fit results agree within statistics with the known inputs, indicating a lack of bias in our signal estimation technique. We have also verified that our method is insensitive to the signalto-background ratio of the sample over a range corresponding to 0.2 − 5.0 times that estimated in the data.   Table VI presents the percentage estimates for statistical and systematic uncertainties from different sources. Figure 9 presents a plot of each uncertainty source as a function of observed multiplicity.

A. Statistical Uncertainties
Statistical uncertainties are returned from the MINUIT package used in our fitting for both data and MC samples. These uncertainties include contributions from the CR background in our fitting procedure, and our procedure includes the CR background systematic uncertainty in a contribution to the total statistical uncertainty. Both data and MC statistics contribute substantially to the overall uncertainties in our data, as shown in Fig. 9.

B. Short Track Efficiency Uncertainties
The dominant systematic uncertainty originates from the differences in the efficiency between data and simulation for reconstructing short-length hadron tracks. The overall efficiencies of the Pandora reconstruction algorithms are a strong function of the number of hits of the tracks, with a plateau not being reached until of order of several hundred hits. The inclusive efficiencies for reconstructing protons or pions at the 15 collection plane hit threshold is estimated to be ε = 0.45 ± 0.05. The absolute efficiency value is not used in this analysis, but we use this estimate to conservatively assign a mean efficiency uncertainty of δ = 15%.
We then estimate the effect of an efficiency uncertainty on multiplicity by the following procedure: Consider a track in an event in a multiplicity bin N. If one lowers the tracking efficiency by the factor 1 − δ , then there is a 1 − δ probability that the track remained reconstructed and the event stayed in that multiplicity bin, and a probability δ that the track would not have been reconstructed and that the event would thus have a lower multiplicity. If the overall multiplicity is N, with N − 1 short tracks and one long track, and each track's reconstruction probability is reduced by a factor 1 − δ , then an overall fraction of events (1 − δ ) N−1 will remain in the bin, and a fraction 1 − (1 − δ ) N−1 will migrate to lower multiplicity bins. The fraction of tracks that migrate to multiplicity N < N from bin N, f (N ; N, δ ), is given by binomial statistics: We use this result to generate the expected observed CPMD in simulation that would emerge from lowering the tracking efficiency by the factor 1 − δ compared to the default simulated CPMD. The difference between the two distributions is then taken as the systematic uncertainty assigned to short track efficiency, with the assumption that the effect of increasing the default efficiency by a factor 1 + δ would produce a symmetric change. Table VII summarizes this study for the three GENIE models used. The observed multiplicity = 1 probability increases because of "feed down" of events from higher multiplicity, due to the lowered efficiency, mainly from observed multiplicity = 2. The other observed multiplicity probabilities decrease accordingly. The largest effects are in high multiplicity bins because the loss of events from lowering the efficiency by the factor (1 − δ ) varies as (1 − δ ) N−1 for multiplicity bin N. Monte Carlo simulations show that "fake tracks" that could move events to higher multiplicity are rare. We have observed no statistically significant differences in the shape distributions after adjusting the efficiency by the constant per-track factor implied by the pull factor.

C. Long Track Efficiency Uncertainties
To first order, the efficiency for reconstructing tracks with length > 75 cm is not expected to affect the observed multi-   plicity distribution, as it is common to all multiplicities and cancels in the ratio when forming observed multiplicity probabilities. At second order, however, a multiplicity dependence that changes the distribution of observed multiplicity without affecting the overall number of events is possible. A plausible model for this is that higher multiplicity in an event helps VII: Relative change in observed multiplicity probabilities corresponding to a −15% uniform reduction in short charged particle tracking efficiencies for three GENIE models: default, MEC, and TEM. The missing entry for observed multiplicity 5 in TEM is due to no event being generated with that observed multiplicity.
Observed multiplicity Pandora better define a vertex, and thus increases the chance that the event passes the ν µ CC selection filter. We estimate the size of this effect by comparing the efficiencies obtained with the Pandora package for simulated quasielastic final states in which both the proton and muon are reconstructed, to charged pion resonance final states in which the proton, pion, and muon are all reconstructed. From this study we conclude that the efficiency for finding the muon in final states where all charged particles are reconstructed could be up to 3% higher for charged pion resonance events (observed multiplicity 3) than quasi-elastic events (observed multiplicity 2). We then assume, for the purpose of uncertainty estimation, that this relative enhancement seen for higher observed multiplicity events in the MC is absent in the data.  VIII: Relative change in observed multiplicity probabilities corresponding to increasing the conditional probability for reconstructing the long track by 3% for each additional track found in the event, as suggested by Pandora studies of QE and charged pion resonance production for three GENIE models: default, MEC, and TEM. The missing entry for observed multiplicity 5 in TEM is due to no event being generated with that observed multiplicity.
Observed multiplicity

D. Background Model Uncertainties
In the signal extraction fitting procedure, two conditional parameters (α ν and α CR ) were extracted from the Monte Carlo simulation and off-beam data. To calculate the systematic uncertainties on these parameters, their values were varied by ±1σ of their statistical uncertainty. Those values were propagated in the observed charged particle multiplicity distribu-tion. We also extracted the α ν and α CR values separately from the GENIE default, GENIE+TEM, and GENIE+MEC models. The effect from this systematic variation were found to be very small.

E. Flux Shape Uncertainties
Variations in flux can be parameterized by where Φ (E ν ) is the neutrino flux at neutrino energy E ν and ∆ (E ν ) is the fractional uncertainty in the flux at that energy. An energy-independent ∆ (E ν ) has no effect on the observed multiplicity distributions as this measurement is independent of absolute normalization. On the other hand, raising the high energy flux relative to the low energy flux could enhance the contributions of higher multiplicity resonance and DIS processes. We confine ourself to considering highly correlated energy-dependent shifts, denoted as ∆ i (E ν ) for i = 1 − 6 via an approximate procedure that should be conservative. These shifts, shown in Fig. 10, are allowed to modify the BNB flux within uncertainties determined by the MiniBooNE collaboration [20]. The first two variations simply shift all flux values up (∆ 1 (E ν )) or down (∆ 2 (E ν )) according to the flux uncertainty envelope. The next two enhance the high energy flux (∆ 3 (E ν )) or low energy flux (∆ 4 (E ν )) linearly with neutrino energy, with the variation taken to be zero at the average neutrino energy. The final two variations enhance high energy flux (∆ 5 (E ν )) or low energy flux (∆ 6 (E ν )) logarithmically with neutrino energy, with the variation taken to be zero at the average energy. As expected, shifts that are positively correlated across all energies produce negligible differences, but even shifts that produce sizable distortions between high and low energies contribute systematic uncertainties that are small.
(GeV)  The measured charge from muon-induced ionization can vary within the detector volume due to the finite probability for drifting electrons to be captured by electronegative contaminants in the liquid argon. This capture probability can be parameterized by an electron lifetime τ. We perform our analysis on simulation with two lifetimes that safely bound those measured during detector operating conditions, τ = 6 msec and τ = ∞ msec. The resulting distribution of percentage uncertainty as a function of multiplicity in Fig. 9 shows that the electron lifetime uncertainties minimally affect the multiplicity.

G. Other Sources of Uncertainty Considered
A systematic comparison was performed on all kinematic quantities entering this analysis between off-beam CR data and the CR events simulated with CORSIKA. No statistically significant discrepancies were observed between event selection pass rates applied to off-beam data and MC simulation.
A check of possible time-dependent detector response systematics was also performed by dividing the data into two samples and performing the analysis separately for each sample. Differences between the two samples are consistent within statistical fluctuations.
The data are not corrected for ν µ NC, ν e ,ν e , orν µ backgrounds. An assumption is made that the Monte Carlo simulation adequately describes these non ν µ CC backgrounds. Section IX A shows that these backgrounds, in total, are expected to be less than 10% of the final sample; their impact on the final distributions is generally small.

A. Observed Charged Particle Multiplicity Distribution
Following the implementation of the signal extraction procedure and verification through closure test on MC events, we execute the same maximum likelihood fit on data. Table IV lists the values of the fit parameters obtained for the data; and Table IX lists the number of neutrino events in different multiplicity bins for the data. While our method does not require this to be the case, we note that the fitted PH and MCS test probabilities P(PH), Q(PH), P(MCS), and Q(MCS) agree in data and simulation within statistical uncertainties. This provides evidence that the simulation correctly describes the muon PH and MCS tests used in the analysis.
Area normalized, bin-by-bin fitted multiplicity distributions from three different GENIE predictions overlaid on data are presented in Fig. 11 where data error bars include statistical uncertainties obtained from the fit and the MC error bands include MC statistical and systematic uncertainties that are listed in Table VI added in quadrature.
In general the three GENIE models agree within uncertainties with one another, and agree qualitatively with the data. There are indications that GENIE overestimates the mean charged particle multiplicity relative to the data. We emphasize that no tuning or fitting has been performed to this or any of the other kinematic distributions.

B. Observed Kinematic Distributions
A key technical feature of our analysis entails performing tests on the pulse height and multiple Coulomb scattering behavior of hits on the long contained track in each event. This allows a categorization of events in each multiplicity into four categories according to whether the long track passes or fails the PH and MCS tests: (PASS, PASS), (PASS, FAIL), (FAIL, PASS), and (FAIL, FAIL). We have shown that the (PASS, PASS) category is "neutrino-enriched" and the (FAIL, FAIL) category is cosmic-ray-dominated. The mixed cases (PASS, FAIL) and (FAIL, PASS) provide samples with intermediate signal-to-background ratios.
Our fit to the distribution of the eight event categories in on-beam and off-beam data allows us to estimate the number of neutrino eventsN νi and the number of corresponding background CR eventsN CRi for each observed multiplicity i. Oncê N νn andN CRn are established, we can obtain a prediction for the content of any bin k of any kinematic quantity X i j associated with track j in an observed multiplicity i event in any (PH, MCS) test combination: Herex ν,i j (PH, MCS) k is an area-normalized histogram of X i j for "true neutrino events" in a given category obtained from a "MC" sample, andx CR,i j (PH, MCS) k is an area-normalized histogram of X i j for CR events obtained from off-beam data. This distribution can be compared to the corresponding one for data in each category, data(X i j , PH, MCS) k .
In short, we assume that the observed distribution of events consists of a mix of neutrino events plus CR events. The proportions of the mix in each category are fixed by the output of our fit, which, by construction, constrains the normalization of the model to equal that of the data. We emphasize that only the PH and MCS tests have been used to extract the neutrino interaction signal sample; no information from any quantity X i j is used.

C. Checks on Distributions lacking Dynamical Significance
Several kinematic properties of neutrino interactions depend only weakly on the neutrino interaction model; these include the reconstructed vertex positions, the initial and final coordinates of the long track, and the azimuthal angles of individual tracks. These distributions provide checks on the overall signal-to-background separation provided by the testcategory fits and flux and detector modeling. They also test for differences between the modeling of neutrino events, which depend on the GEANT detector simulation, and CR events, which use the off-beam data and thus do not depend on detector simulation.
As an example, we show the observed distributions for the selected vertex y position for the candidate muon track from the full selected sample in Fig. 12. For this and all subsequent distributions, the on-beam data events are indicated by plotted points with statistical error bars. The model prediction is shown by a colored band (red for GENIE default, green for GENIE+TEM, and blue for GENIE+MEC) with the width of the band indicating the correlated statistical plus efficiency systematic uncertainty from using common N ν,n , N CR,n values for all bins of all distributions of a given multiplicity bin n. The CR contribution to a distribution in a given category is shown by the shaded cyan region. For example, Fig. 12 com-pares the on-beam data to GENIE default MC sample and also shows the CR background.
The signal-enriched (PASS, PASS) category for vertex y has the nearly flat distribution expected for a neutrino event sample with a small CR background. Note that in our selection, we only allow candidate muon tracks initial y position < 70 cm. This cut rejects many cosmic rays that produce a downward trajectory in the final selected sample. The remaining background is dominated by cosmic rays with an apparent upward trajectory. This can be seen in the background-enriched sample (FAIL, FAIL) in the vertex y distribution where a peak at negative y values corresponds to "upwards-going"CR. Figure 13 shows the distribution of azimuthal angle φ , defined in the plane perpendicular to the beam direction, of the muon candidate track for the full selected sample. The CRdominated (FAIL, FAIL) category shows the expected peaking at φ = ±π/2 from the mainly vertically-oriented CR. The asymmetry in the peak's structure is due to the requirement on vertex y position described previously in Sec. V D. By contrast the signal-enriched (PASS, PASS) category has the nearly flat distribution expected for a neutrino event sample with a small CR background.
Similar levels of agreement exist between data and simulation for distributions of the event vertex x and z positions, for the (x, y, z) position of the end point of the muon track candidate, and for the azimuthal angles of individual tracks in multiplicity 2 and 3 topologies. We thus conclude that the simulation and reconstruction chain augmented by our method for estimated CR backgrounds satisfactorily describes features of the data that have no dependence on the neutrino interaction model.

D. Dynamically Significant Distributions
Events with N reconstructed tracks have potentially 4N dynamically significant variables−the components of each parti-  cle 4-vector−which will have distributions that depend on the neutrino interaction model. Azimuthal symmetry of the beam eliminates one of these, leaving 3, 7, 11, and 15 dynamically significant variables for multiplicities 1 − 4, respectively. In the following, we use the notation X i j to label a dynamical variable x associated with track j in an observed multiplicity i event. For example, cos θ 11 describes the cosine of polar angle distribution of the only track in multiplicity 1 events, while L 22 would describe the length of the second (short) track in multiplicity 2 events. The notation with three subscripts, X i jk , represents a distribution of the difference in variable x associated with tracks j and k in an observed multiplicity i event.
For one-track events, three variables exist. We use the observed length L 11 as a proxy for kinetic energy, and the cosine of the scattering angle with respect to the neutrino beam direction cos θ 11 . The azimuthal angle φ 11 has no dynamical significance and must be uniformly distributed due to the cylindrical symmetry of the neutrino beam.
Since the particle mass is not determined in our analysis, we are free to introduce a third dynamically significant quantity that is sensitive to particle mass, which we take to be whereŝ 11 is a unit vector parallel to the track direction at the event vertex, andt 11 is a unit vector that points from the start of the track to the end of the track in the detector. The variable Θ i j measures the angular deflection of a track over its length due to multiple Coulomb scattering. Its dependence on track momentum and energy differs from that of track length. For most of the MicroBooNE kinematic range, we expect light particles (π and µ) to scatter more, and thus produce a broader sin Θ i j distribution, than protons over the same track length. Figure 14 shows the distributions of L 11 , cos θ 11 , and sin Θ 11 , from the neutrino-enriched sample compared to the GENIE default model. Figure 15 presents the L 11 distribution for the GENIE+MEC and GENIE+TEM models. This is the distribution where the agreement between data and GE-NIE+TEM model, compared to the agreement between data and the other two models, is largest. Figure 16 presents the cos θ 11 distribution for GENIE+MEC and GENIE+TEM models. This is the distribution where the agreement between data and the GENIE default compares least favorably than to the GENIE+MEC and GENIE+TEM models.
For brevity in the following, except where noted, we only show comparisons of data to predictions from the GENIE default model.
Comparisons to GENIE+TEM and GE-NIE+MEC show qualitatively similar levels of agreement. Differences for specific distributions can be examined in terms of the χ 2 test statistic values in Table X. For two track events, seven dynamic variables exist. These include properties of the long track that parallel the choices for one-track events, L 21 , cos θ 21 , and sin Θ 21 , similar quantities for the second track, L 22 , cos θ 22 , and sin Θ 22 , plus a quantity that describes the correlation between the two tracks in the event, which we take to be the difference in azimuthal variables φ 221 = φ 22 − φ 21 . Since track 2 can exit the detector, the meaning of L 22 and sin Θ 22 differ somewhat from L 21 and sin Θ 21 . Two other two-track correlated variables of interest, which are not independent, are the cosine of the opening angle, cos Ω 221 = cos θ 21 cos θ 22 + sin θ 21 sin θ 22 cos (φ 22 − φ 21 ) , (31) and the cosine of the acoplanarity angle withẑ a unit vector in the neutrino beam direction andŝ 21 is a unit vector parallel to the first track direction at the event vertex. For the scattering of two initial state particles into two final state particles (2 → 2), one expects from momentum conservation φ 221 = ±π and cos θ A = 0. Deviations of φ 221 from ±π or of cos θ A from 0 could be caused by undetected tracks in the final state, from NC events in the sample, or from effects of final state interactions in CC events. The opening angle serves a useful role in identifying spurious two-track events that result from the tracking algorithm "breaking" a single track into two tracks, most commonly in cosmic ray events. Broken tracks produce values of cos Ω 221 very close to −1. Figures 17 and 18 show the distributions of (L 21 and L 22 ) and (cos θ 21 and cos θ 22 ) from the neutrino-enriched sample, compared to the GENIE default model. Figure 19 presents the distributions of cos θ 21 using GENIE+MEC and GENIE+TEM models. This is the distribution where the agreement between data and GENIE default model, compared to the agreement between data and the other two models, is largest. Figures 20 and 21 show the distributions of (sin Θ 21 and sin Θ 22 ) and (φ 22 −φ 21 and cos Ω 221 ) from the neutrino-enriched sample, compared to the GENIE default model. We have performed a test where we remove all events in the first bin of Fig. 21; we see no changes in the level of agreement between data or model in other kinematic distribution comparisons, and no statistically significant shifts in the observed multiplicity distributions.

E. χ 2 Tests for Kinematic Distributions
We quantify agreement between model and observation through use of χ 2 tests on the kinematic distributions described in Sec. VIII D. Ensemble tests have established the validity of the use of the χ 2 criterion. We use only the "neutrino-enriched" sample of events in which the candidate muon passes both the PH and MCS tests. Data are binned     into histograms, with a bin k for a variable x i j , d i jk , and compared to model predictions constructed by assuming that the number of events in a bin k of a variable X i j , m i jk , consists of contributions from neutrino and CR background contributions. We shorten the notation in Eq. 29 to where M ν,i and M CR,i are the number of neutrino and CR events, respectively, predicted to be in the neutrino-enriched category for multiplicity i (as described in Sec. VIII B); and x ν,i jk andx CR,i jk the fraction of neutrino and CR events, respectively, falling in the k bin for variable x i j as predicted by the GENIE model and the off-beam CR sample, respectively. Thex ν,i jk andx CR,i jk are shape distributions normalized to one: We then construct a χ 2 for x i j using a Poisson form appropriate for the low statistics in many bins:    Table X summarizes the results of these χ 2 comparison tests for 21 independent kinematic variables to the three GENIE models. Only bins with at least one data event and one model event were used in the calculation of χ 2 . The number of degrees of freedom associated with the χ 2 test was set equal to the number of bins used for that histogram minus one to account for the overall normalization adjustment. We note here that these tests for consistency are defined at the level of statistical uncertainties only; systematic uncertainties are not incorporated into the χ 2 terms.
We summarize here salient features of Table X as follows: All three models consistently describe the data, with summed χ 2 per degree-of-freedom (χ 2 /DOF) of 228.1/216, 216.9/216, and 229.6/216, respectively, and corresponding p-values of P χ 2 = 27%, 47%, and 25% for GENIE default, MEC, and TEM, respectively. The total χ 2 after including all dynamic and non-dynamic variable distributions is 714/652. No tune of GENIE is superior to any other with any meaningful statistical significance for the distributions we have considered. The acceptable values of χ 2 are consistent with the hypothesis that the combination of a GENIE event generator, the MicroBooNE BNB flux model, and the MicroBooNE GEANT-based detector simulation satisfactorily describe the properties of neutrino events examined in this analysis in a shape comparison. All elements of the MicroBooNE analysis chain thus appear to be performing satisfactorily; and no evidence exists for missing systematic effects that would produce data-model discrepancies outside the present level of statistics.
Aggregating Table X different ways uncovers no significant discrepancies. The χ 2 tests on leading track cos θ and sin Θ yield satisfactory results for all multiplicities. Combined χ 2 /DOF for all distributions associated with a particular multiplicity likewise exhibit adequate agreement. The most poorly described single distribution is that for the length of the muon candidate in multiplicity 1 events. The P χ 2 , while acceptable, are 16% and 8% for the GENIE default and GENIE+MEC, respectively. The GENIE+TEM model has P χ 2 = 53%.
The χ 2 values for different distributions in a given multiplicity are calculated using the same events, which gives rise to concerns about correlations between different distributions. We have performed studies that verify that the χ 2 values would be highly correlated if the model and data disagreed by an overall normalization, but that otherwise the χ 2 tests on different distributions exhibit independent behavior, even when the same events are used. The P χ 2 values for different distributions do not cluster near 0 or 1, which is consistent with the view that the projections display approximately independent statistical behavior.
In summary, all GENIE models successfully describe, through χ 2 tests, the shapes of a complete set of dynamically significant kinematic variables for observed charged particle multiplicity distributions 1, 2, and 3. The statistical powerthe highest precision afforded by the available statistics with which the predictions can be tested-of these tests from the overall data statistics available corresponds to approximately 4%, 7%, and 20% for multiplicity 1, 2, and 3, respectively.

F. χ 2 Tests for Multiplicity Distribution
While we find satisfactory agreement between GENIE models and kinematic distribution shapes using χ 2 tests that in-corporate only statistical uncertainties, the situation differs for the overall multiplicity distribution. Here, we find statistical χ 2 M /DOF= 30/4, 22/4, and 28/4 for the default, MEC, and TEM GENIE models, respectively. However, in the case of multiplicity, a significant systematic uncertainty exists for tracking efficiency that must be taken into account before any conclusion can be drawn.
We incorporate a tracking efficiency contribution to the χ 2 test by defining Here D i is the number of neutrino events estimated by the signal extraction procedure (Sec. VIII B), and σ i is the estimated uncertainty on D i using the signal extraction procedure. For multiplicity 3 and higher, the uncertainty on D i is purely statistical as the CR background becomes negligible. The quantities M i are the number of events in the MC sample with multiplicity i. Finite statistics in the MC sample are incorporated by interpreting the M i as Poisson fluctuations about their true valuesM (0) i in the third term of Eq. 38. This analysis does not absolutely normalize MC to data, hence the relative normalization of data to MC is allowed to float via the parameter K in the first term of Eq. 38. The normalization constant K, while not used directly in the model test, is consistent with the predicted value from the default GENIE model.
As discussed in Sec. VII B, changing the per-track efficiency by a constant fraction δ in the model would shift events between multiplicities according tô For the nominal model used in the MC simulation δ = 0. As discussed in Sec. VII B we estimate the uncertainty on δ to be 15%, and we introduce this into χ 2 M through the "pull term" (δ /0.15) 2 .
We minimize χ 2 M with respect to the tracking efficiency pull parameter δ , the MC-to-data normalization K, and the five MC We find that a satisfactory χ 2 value can be obtained for the multiplicity distribution itself, albeit at the cost of a ≈ 2σ pull in the parameter δ .
The following list summarizes qualitative expectations for components of observed multiplicities from particular processes. These components can include contributions from the primary neutrino-nucleon scatter within the nucleus and secondary interactions of primary hadrons with the remnant nucleus. Secondary charged particles are usually protons, which are expected to be produced with kinetic energies that are usually too low for track reconstruction in this analysis. However, more energetic forward-produced protons from the upper "tail" of this secondary kinetic energy distribution may contribute. Note that the particle-type-dependent kinetic energy thresholds for charged particles entering our sample range from 31 MeV for a π ± to 69 MeV for a proton.
• Multiplicity > 3, mainly predicted to be "DIS events" in which at least three short tracks are reconstructed. "DIS" is the usual term for multi-particle final states not identified with any particular resonance formation.
• Multiplicity = 3, mainly predicted to be µ − pπ + events from ∆ resonance production in which all three tracks are reconstructed. "Feed down" from higher multiplicity would be small due to the relatively small DIS cross section at MicroBooNE energies.
• Multiplicity = 2, mainly predicted to be QE µ − p events and resonant µ − pπ 0 events in which the proton is reconstructed, with a sub-leading contribution from "feed down" of resonant charged pion production events where one track fails to be reconstructed.
• Multiplicity = 1, mainly predicted to be "feed down" from QE µ − p and µ − pπ 0 events in which the proton is not reconstructed, with contributions from other higher multiplicity topologies in which more than one track fails to be reconstructed. Figure 28 illustrates these expectations from GENIE. We note that, as expected, the three-track topology is dominated by resonant pion production in the default GENIE model, while the two-track and one-track topologies are QEdominated with non-negligible resonance feed-down. The coherent pion production process (ν µ + Ar → µ − π + + Ar) denoted by "CCCohP" in this figure, as well as NC and ν e and ν e scattering, only lead to small contributions.
Our observation of discrepancy of data compared to simulation in three-track compared to two-track topologies, shown in Fig. 11, is qualitatively similar to the low ν µ CC pion cross sections compared to GENIE reported by MINERvA [38] using hydrocarbon targets at the somewhat higher neutrino energy from the Fermilab Neutrinos at the Main Injector (NuMI) beam. The T2K experiment [39] also reports a low pion production cross section relative to GENIE expectations using water targets in a neutrino beam with comparable energy to the BNB. However, MiniBooNE measured a charged pion production rate more in agreement with GENIE using mineral oil as a target in the same Fermilab BNB as used by Micro-BooNE [40].
MicroBooNE also observes more one-track events than GE-NIE predicts, as shown in Fig. 11. This corroborates Ar-goNeuT's observation that approximately 35% of neutrino interactions on argon targets with no pions detected in the final state also contained no observable proton [41].
While our observed multiplicity distribution disagrees with GENIE expectations and shows consistency with a number of other experiments, we cannot, as noted in Sec. VIII F definitively exclude an alternate explanation of the discrepancy in terms of a tracking efficiency error at this time.
Our kinetic energy thresholds limit acceptance in such a way that protons produced in FSI may not significantly contribute to the observed CPMD. Furthermore, our analysis requires a forward-going long contained track as a muon candidate, which restricts the final state phase space. Also, our analysis makes use of fully automated reconstruction. Therefore, results of this analysis should not be directly compared to the low energy proton multiplicity measurement reported by ArgoNeuT [17].

B. GENIE Predictions for Kinematic Distributions
Kinematic distributions for fixed multiplicity suffer much less from tracking-related systematic uncertainties than the multiplicity probabilities; hence GENIE expectations for the shapes of kinematic distributions can be compared directly to data. The MEC and TEM tunes of GENIE primarily change the normalization of QE-like event topologies relative to resonance type topologies, and secondarily modify properties of low energy final state protons that would usually not satisfy our acceptance criteria. Shape comparisons would thus not be expected to differentiate MEC and TEM from the GENIE default, and we have verified this expectation with our χ 2 tests. Accordingly we confine the following discussion to the default GENIE tune. Figure 29 shows the predictions for reconstructed L 11 , cos θ 11 , and sin Θ 11 from the neutrino-enriched sample, using the GENIE default model. The muon track candidate is only mildly affected by the details of the recoiling hadronic system, and thus QE, RES, and DIS production produce similar shape contributions to L 11 , cos θ 11 , and sin Θ 11 . Figures 30, 31, 32, and 33 present the distributions of reconstructed (L 21 and L 22 ), (cos θ 21 and cos θ 22 ), (sin Θ 21 and sin Θ 22 ), and (φ 22 − φ 21 and cos Ω 221 ) respectively from the neutrino-enriched sample, using the GENIE default model. There is again minimal difference between QE, resonance, and DIS channels in track length, cos θ , or sin Θ for the leading track. However, the QE channel produces contributions to the distributions in cos θ 21 and cos θ 22 that are considerably less forward-peaked than the resonance channel contributions. Distributions of these quantities in the data appear to be consistent with this picture. We also note that the sin Θ 22 distribution receives a contribution from QE scattering peaked at small values, consistent with expectations for a proton, and a broader distribution more similar to that of the leading muon track candidate that is consistent with the hypothesis that a charged pion can be reconstructed as the second track in resonance contributions.

X. SUMMARY
We have completed an analysis that compares observed charged-particle multiplicities and observed kinematic distributions of charged particles for fixed multiplicities in a restricted final state phase space for neutrino scattering events in argon to predictions from three GENIE tunes processed through the MicroBooNE simulation and reconstruction chain. Our analysis takes into account statistical uncertainties in a rigorous manner, and estimates the impact of the largest expected systematic uncertainties. We observe that all elements of the MicroBooNE measurement chain−detector performance, data acquisition, event reconstruction, Monte Carlo event generator, detector simulation, and flux modeling−perform well.
With particle-type-dependent kinetic energy thresholds of 31 MeV for π ± and 69 MeV for protons, we find all three GENIE tunes consistently describe data in the shapes of 26 different kinematic distributions at fixed multiplicities. GE-NIE appears to over-predict the number of three-track events in data that would be expected from resonant pion production, and to under-predict the number of one-track events; however, we cannot rule out a higher than expected tracking efficiency uncertainty as an alternative explanation for these observations. Our study thus empirically supports the use of GENIE in describing single-process (quasi-elastic, resonance) neutrino scattering on argon, but not the predictions for the relative contributions of different processes to the overall cross section. We find no significant differences at this stage in the experiment between the default GENIE tune or tunes that add MEC or TEM. Use of any of the three GENIE tunes for future MicroBooNE analyses, or for physics studies of inclusive final states performed for the SBN and DUNE experiments, receives empirical validation from this work.
As part of this analysis, we have developed a data-driven cosmic ray background estimation method based on the energy loss profile and multiple Coulomb scattering behavior of muons. Within the available Monte Carlo statistics, we have shown that this method provides an unbiased estimate of the number of neutrino events in a pre-filtered sample, and, given current statistical precision, it is independent of the signal-tobackground level, final state charged particle multiplicity, and other kinematic properties of the final state particles. This method can be applied to a broad range of charged current process measurements.
Significant improvements to MicroBooNE neutrino interaction property measurements are anticipated in the future through incorporation of nearly an order-of-magnitude more statistics, more fully developed reconstruction tools−including momentum reconstruction, particle identification, and lower kinetic energy thresholds for tracking−and the availability of a recently installed external cosmic-ray tagger to the detector. This material is based upon work supported by the following: the U.S. Department of Energy, Office of Science, Offices of High Energy Physics and Nuclear Physics; the U.S. National Science Foundation; the Swiss National Science Foundation; the Science and Technology Facilities Council of the United Kingdom; and The Royal Society (United Kingdom). Additional support for the laser calibration system and cosmic ray tagger was provided by the Albert Einstein Center for Fundamental Physics. Fermilab is operated by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the United States Department of Energy.