The extendedtrack reconstruction for MiniBooNE
Abstract
The Booster Neutrino Experiment (MiniBooNE) searches for oscillations using the neutrino beam produced by the FNAL Booster synchrotron. The array of photomultiplier tubes (PMTs) lining the MiniBooNE detector records Cherenkov and scintillation photons from the charged particles produced in neutrino interactions. We describe a maximum likelihood fitting algorithm used to reconstruct the basic properties (position, direction, energy) of these particles from the charges and times measured by the PMTs. The likelihoods returned from fitting an event to different particle hypotheses are used to categorize it as a signal event or as one of the background processes, in particular charged current quasielastic scattering and neutral current production. The reconstruction and event selection techniques described here can be applied to current and future neutrino experiments using similar Cherenkovbased detection.
keywords:
Pacs:
,
,
,
,
,
1 Introduction
Cherenkov detectors of have played a key role in the establishment of the phenomenon of neutrino oscillations [?, ?, ?]. These detectors typically consist of a large volume of a homogeneous transparent medium (water or mineral oil) with a high index of refraction surrounded by an array of photomultiplier tubes (PMTs). Cherenkov photons produced by charged particles emerging from the neutrino interactions and traversing the medium with (where is the index of refraction and ) are detected by the PMTs. The photons are emitted at an angle relative to the track direction, where . The radiation is azimuthally symmetric about the track direction, resulting in a ringlike pattern that can be identified on the PMT array. The quantity, spatial distribution, and arrival times of these photons provide information on the location, direction, and energy of the particle. We refer to the extraction of such information from the charge and time measurements of the PMTs as “reconstruction.” An analysis of the ring profile can also provide information on the identity of the particle and, if multiple rings are detected, the number of particles.
In this paper, we discuss the event reconstruction algorithms used in the Booster Neutrino Experiment (MiniBooNE), which searches for an excess of interactions indicative of oscillations using the Fermilab neutrino beam. The MiniBooNE detector is a sphere of radius 610.6 cm filled with Marcol 7 mineral oil () and divided into two optically isolated regions by an opaque shell of radius 575 cm, concentric with the sphere. The surface of the inner “main” region is instrumented with an array of 1280 inwardfacing 8 in. PMTs which detect the light produced by the neutrino interactions. The outer “veto” region is instrumented with 240 PMTs and is used to tag charged particles that enter the detector from outside (e.g., cosmic muons) or that exit the main region. Though no scintillator was added to the Marcol 7, it scintillates weakly, resulting in the production of delayed isotropic light for particles with sub and superCherenkov velocities. A more detailed description of the experiment can be found in Refs. [12] and [13].
2 Neutrino Interactions at MiniBooNE
The neutrino beam at MiniBooNE results in interactions with relatively low outgoing particle multiplicity. The largest interaction channel is the charged current quasielastic (CCQE) process:
(1) 
which accounts for of all neutrino interactions in the MiniBooNE detector. Since the recoil proton of Eq. (1) is typically below Cherenkov threshold, only the outgoing lepton produces significant light. In modeling such events within a reconstruction algorithm, one can therefore consider only the Cherenkov and scintillation light produced by the outgoing lepton. While the recoiling nucleon can produce scintillation light, this additional source of light is not considered in the reconstruction.
A muon in the MiniBooNE detector exhibits minimumionizing behavior through most of its path with little chance of radiative energy loss. In contrast, electrons typically induce electromagnetic showers, resulting in additional electrons and positrons that emit Cherenkov light. Thus, muons and electrons have significantly different Cherenkov ring patterns. This difference is the basis for distinguishing a muon track from an electron track, and the reconstruction algorithm has a model for each.
The next most abundant process in MiniBooNE is single pion production, which occurs primarily via resonance or coherent scattering:
(2) 
where denotes a nucleon in the case of resonance production and a nucleus in the case of coherent production. Of particular interest is neutral current (NC) production. The two photons from decay induce electromagnetic showers indistinguishable in MiniBooNE from those induced by electrons. Since the recoil nucleon/nucleus is typically below Cherenkov threshold, these NC events are welldescribed in the reconstruction algorithm by two electronlike tracks pointing back to a common vertex. The presence of two distinct Cherenkov rings allows for the separation of events from electron ( CCQE) events. However, misidentification (for example, due to a large energy asymmetry between the two photons or due to a small photon opening angle which leaves the Cherenkov rings overlapping) accounts for most of the reducible background in the oscillation search. Thus, successful reconstruction and identification of events is critical.
While other channels can be accommodated by the reconstruction algorithm, the appearance oscillation analysis uses only four event models: single electron track, single muon track, two tracks, and two tracks with a invariant mass.
3 The Detector Response
The Marcol 7 mineral oil used in MiniBooNE exhibits a rich array of optical phenomena despite its m extinction length near the peak of PMT sensitivity ( nm). Cherenkov and scintillation light production is accompanied by photon absorption, fluorescence (with several possible excitation/emission spectra and lifetimes), and Rayleigh and Raman scattering. The left plot in Figure 1 summarizes the rates of various processes as a function of wavelength. The cumulative extinction rate is shown as the black line. In the nearultraviolet region (), fluorescence processes dominate the extinction rate. The fluorescence lifetimes range from . In the visible region (), the dominant processes are Rayleigh scattering and absorption.
The 1520 8 in. PMTs in MiniBooNE are of two types: 322 model R5912 PMTs and 1198 model R1408 PMTs, both from Hamamatsu. All of the R5912 PMTs are located in the main PMT array. The PMT time response (particularly the latepulsing behavior) and the variation of PMT efficiency with incident angle have been characterized in external measurements using a pulsed LED [14]. The PMTs in situ have been studied using a laser calibration system that produces subnanosecond bursts of nearly isotropic light in the detector.
The right panel of Figure 1 shows the times of PMT hits recorded in laser calibration events relative to their expected times. The colored histograms come from Monte Carlo simulation, with the red curve showing PMT and electronics timing effects (including pre and latepulsing), the green curve including photon scattering, and the blue curve adding photon reflections from the PMTs and the surface of the main detector region. The resulting time structure matches well with that seen in data (black points). We note that only the earliest photoelectron’s time is reported when multiple photoelectrons are present in a hit. (A fresh hit can be seen after an electronics dead time of ns.) The reported charge, however, does account for additional photoelectrons.
Absorption, scattering, reflections, and fluorescence processes influence both the topology and the time structure of hits. Additionally, scintillation light, emitted with a lifetime of , adds intrinsically delayed light to the prompt Cherenkov photons. The event models within the reconstruction algorithm account for all these phenomena. In the discussion, we refer to Cherenkov and scintillation light that travels undisturbed from the particle track to the PMT as “direct” light. Light that undergoes fluorescence, scattering, or reflection we term “indirect” light.
A detailed GEANT3based [15] Monte Carlo simulation of the MiniBooNE detector serves as the central tool for developing the reconstruction algorithm’s predictive models. The accuracy of the models influences performance through resolution and biases in extracted track parameters and through the particle separation capabilities of the likelihoodbased hypothesis testing. To balance the need for likelihoods that are as accurate as possible with the desire to contain the computational demands, detector and track model information extracted from the Monte Carlo simulation is tabulated, approximated, or parametrized wherever possible while minimizing any sacrifice in performance. Also, thanks to the excellent stability of the detector’s optical properties, a single configuration of the reconstruction is sufficient for use on all MiniBooNE physics data.
4 The Single Track Model
The reconstruction of a single track in the detector, whether an electron or muon, is based on a model with seven parameters:

the starting point:

the starting time:

the direction: ,

the kinetic energy: .
The starting point and time are referred to as the “vertex” of the event. The direction is described in terms of the polar angle defined with respect to the direction of the neutrino beam and the azimuthal angle about this beam axis. For the electron model, the stochastic variations in the electromagnetic showers are modeled in average – no attempt is made to account for eventbyevent variations in the shower profile. This is also true of the much smaller variations in the muon propagation (e.g., straggling). This simplification allows the track to be fully specified from its initial properties summarized in the seven parameters above. We refer to this vector of parameters as .
The observables of an event are the measurements from the 1280 PMTs in the main region of the detector. For each PMT, we have:

a bit indicating whether the tube registered a hit

if the tube was hit, the measured charge of the hit

if the tube was hit, the measured time of the hit.
Assuming that the PMTs behave independently, the likelihood for an observed set of PMT measurements given track parameters can be expressed as:
(3) 
where the products are taken over all unhit and hit PMTs and:

is the probability that the th tube is hit, given parameters .

, are the measured charge and time on the th PMT.

is the probability distribution function (PDF) for the measured charge on the th PMT, given , evaluated at .

is the PDF for the measured time given , evaluated at .
It is convenient to work with the negative logarithm of , and since the charge and timerelated portions of the likelihood are distinct, it is helpful to define:
(4) 
where
(5) 
For brevity, we refer to and as the charge and time likelihoods, respectively, although they are the negative logarithms of the likelihoods.
4.1 The Charge Likelihood
If the number of photoelectrons produced in a PMT is known, then for the observed charge is fully specified in terms of the response properties of the PMT (assumed fixed) without any reference to any other property of the detector. Further, is a Poisson variable whose mean is a function of the track parameters . These considerations motivate a twostep approach to calculating the charge likelihood.
In the first step, using the known optical photon and particle propagation properties of the detector, one determines for a given particle type (electron/muon) and set of track parameters the average number of photoelectrons that a particular PMT should observe. This quantity is referred to as , the “predicted charge” for the th PMT. In calculating the predicted charge, one considers the quantity and angular distribution of light produced by the track, the absorption, scattering, and fluorescence of light as it traverses the mineral oil, the acceptance of the PMT, and anything else that influences the mapping .
For the second step, given in Eq. (5) can be rewritten in terms of the predicted mean charges :
(6) 
where and are now functions of the , with the mapping from the first step left implicit. These two functions depend only on the properties of the PMT response and are, in particular, independent of the happenings within the detector volume. This stepwise approach decouples the track and optical photon model from the PMT response model. The latter is established once using laser calibration data. The laser system provides a source of approximately isotropic and prompt light with controllable intensity. The value of for a particular laser intensity can be determined from the occupancy of the PMTs. The function for each PMT type can be determined by obtaining the distribution of measured charges seen at various intensities . These distributions incorporate all processes associated with the charge response, such as the cascade of charge down the dynode chain, the amplification and digitization of the anode signal, and the conversion of the digital readouts into calibrated charge values. Examples of obtained in this way are shown in Figure 2.
4.2 Calculating the Predicted Charge
We begin our description of the predicted charge calculation with a simple scenario involving a fictitious point source of isotropic direct light. This is then generalized to a line source of isotropic light, appropriate for the scintillation production of a track. We then consider a line source in which the emission has a nontrivial (but azimuthally symmetric) angular distribution, representative of Cherenkov radiation. Finally, we account for the contribution of indirect light (i.e., scattering, reflections, etc.) from both the scintillation and the Cherenkov emission. In order to simplify notation and without loss of generality, the discussion involves fixed but arbitrary PMT location and track parameters .
One needs only two parameters to describe the geometric relationship between an isotropic point source and a PMT. We choose (1) the distance between the source and the PMT and (2) the angle of incidence at the PMT, with corresponding to incidence parallel to the PMT axis. The predicted charge on the PMT can be written
(7) 
where is an eventenergydependent scintillation light yield (i.e., ), is a dependent solid angle factor, is the transmission of the oil and PMT glass as a function of light propagation distance, and is the acceptance of the PMT as a function of the angle of incidence. The transmission depends on the wavelength spectrum of the light source and, thus, is specific in this case to scintillation. The wavelengthdependent PMT photocathode efficiency is also included in . The solid angle factor and PMT angular acceptance are purely geometric, independent of the particulars of the light production.
To generalize to an extended source of scintillation light, as depicted in Figure 3(left), the latter three factors in Eq. (7) become functions of , the distance along the track from its origin. We must also include an emission profile describing the distribution of scintillation production along the track. This profile satisfies . The predicted charge can now be obtained by integrating along the track’s axis:
(8) 
The dependences of , , and on and have been recast as dependences on , as the relationship is fixed for given track parameters and PMT location. Figures 5 and 5 show for 300 muons and electrons, obtained via Monte Carlo simulation. The distribution is relatively flat for muons until the end of the track, where the ionization rate per unit track length becomes larger. (Saturation effects are taken into account in the simulation.) For electrons, the distribution reflects the showering behavior of the track, with rising and falling with the number of shower particles.
We now consider an extended track emitting light with a nontrivial angular distribution, as is the case for Cherenkov radiation. We introduce , the angle to the PMT from the track, as shown in Figure 3(right), and express the predicted charge due to Cherenkov radiation as:
(9) 
This expression differs from its scintillation counterpart by the presence of an angular emission profile . Note that depends on in two distinct ways: the angular profile of the light changes as the track propagates and loses energy, and the angle to the PMT changes depending on which part of the track we are considering.
Figures 7 and 7 show and for simulated 300 muons and electrons. As the muon propagates, loses energy, and approaches the Cherenkov threshold, the rate of Cherenkov radiation per unit track length decreases and the Cherenkov angle becomes smaller. The scattering of the muon, which causes deviations of the track from its original direction, has been included in the Monte Carlo simulation and results in the spread of the angular distribution about the nominal . For electrons, follows the shower profile, like . The presence of the shower particles is readily apparent in the distribution, which becomes substantially wider as the “track” propagates.
In writing Eq. (9), a few simplifications are made. A more rigorous treatment would involve , the differential probability of sending an emitted Cherenkov photon in the direction , which would be integrated over the solid angle subtended by the PMT. By assuming that is constant over the PMT face, the integration can be effected by simply multiplying the differential probability by the solid angle . Also, the azimuthal symmetry of the emission reduces the twodimensional PDF to a onedimensional expression, namely . We take to satisfy
(10) 
for all values of , with the result that a factor of is absorbed into the definition .
4.3 Indirect Light (scattering, fluorescence, reflections)
The above formalism determines the predicted charge for light arriving at the PMTs directly from the track without any redirection. However, the detector has sources of indirect light from scattering, fluorescence, and reflection as discussed in Section 3.
The geometry for indirect light given scintillation (isotropic) emission is shown in Figure 8(left), where an infinitesimal element along the track is situated at radius from the center of the detector and at angle relative to the position of the PMT. The direct light from this track element is simply the integrand of Eq. (8):
(11) 
An analytic expression for the indirect light would involve an elaborate integral over emission angles and scattering points throughout the tank. Rather than attempt this, we observe that the value of such an integral must be proportional to the source strength and must otherwise depend only on the topological variables and . The source strength dependence can be eliminated by forming a ratio of the indirect and direct light predictions:
(12) 
, which we refer to as the scattering table and which we build via the detector simulation, is a property only of the detector optics and the of the track element. With this table, the indirect light contribution can be immediately incorporated into the expression for predicted charge:
(13) 
where the dependence of and on has been made explicit.
For Cherenkov light, the situation is more complex since the light emission is anisotropic. Two additional variables are needed to specify the direction of a vector source relative to the PMT position and the tank center. We make the following nonunique choice for these two variables, as shown in Figure 8(right):

: the angle between (1) the source direction vector and (2) the sourcetoPMT ray (the same as defined elsewhere).

: The angle between (1) the plane containing the tank center, the PMT, and the source, and (2) the plane containing the track and the tank center.
The intensity of the indirect Cherenkov light is normalized with respect to a fictitious isotropic Cherenkov source with predicted charge , so that:
(14) 
The total Cherenkov contribution to the mean predicted charge, including both direct and indirect components is
(15) 
where is given by Eq. (9) and is given by
(16) 
Implicit in this expression is the assumption that the Cherenkov angular emission profile does not change as the track propagates (and loses energy). This simplification has negligible impact, as the processes that result in indirect light tend to destroy the initial emission pattern anyway.
We now have a complete charge prediction, where the total predicted charge for the PMT is given by:
(17) 
4.4 Computation of Integrals
The predicted charge integrals above would be impractical to evaluate with sufficient spatial granularity within a maximum likelihood fit, where they must be evaluated multiple times for every PMT in every event. To mitigate this, all integrations are performed and tabulated beforehand as follows.
The integrand of Eq. (8), which pertains to direct scintillation light, is the product of a “source” factor and an “acceptance” factor
(18) 
As in Eq. (8), the dependence on and is recast through . If varies gradually enough, it can be approximated with a parabolic form:
(19) 
An example of this is shown in Figure 9. The predicted charge then becomes
(20) 
The first integral is identically unity. The remaining two depend on the energy via , but they depend on no other track parameters. This allows one to tabulate the integrals beforehand as a function of , eliminating their evaluation from the actual likelihood calculation.
The parabolic coefficients are obtained by evaluating Eq. (18) at three points along the track: the start of the track (), the midpoint of the Cherenkov emission profile () and twice this distance from the start (). The choice of where to evaluate is somewhat arbitrary; one could choose any three points that sample a large fraction of the range over which light is produced. The extension to indirect isotropic light is straightforward with a redefinition of that incorporates the scattering table:
(21) 
An analogous threeintegral expression can be formed for the Cherenkov predicted charge:
(22) 
As in the scintillation case, the integrals depend on the energy of the track. However, they also depend on two parameters which define the PMTtrack geometry. We choose these parameters to be the vertextoPMT distance, ; and the cosine of the angletoPMT viewed from the vertex, . Labeling the integrals , we obtain
(23) 
Although the integrals depend on three parameters, it is still feasible to tabulate their values ahead of time, eliminating their explicit evaluation within the likelihood calculation.
For indirect Cherenkov light, recall that the scattering table is normalized to a fictitious isotropic Cherenkov source. Thus, the parabolic method used for indirect scintillation can be applied here to give
(24) 
with two energydependent integrals to be tabulated.
4.5 The Time Likelihood
The time portion of the likelihood involves the PDF for the measured PMT hit time given track parameters . Some of the dependence on the track parameters can be eliminated by using the “corrected time” :
(25) 
The expression removes from three terms, namely the starting time of the track , the expected time for light to propagate from the track midpoint to the PMT (with the speed of light in mineral oil given by ), and the time for the particle to propagate from its starting point to the midpoint. The dependence of on the track energy is made explicit. Defining with respect to the track midpoint (as opposed to, say, the start of the track) improves the validity of the simplifications we make below.
Given the substitution, we seek the PDF for the corrected time given arbitrary PMTtrack configurations. The space of configurations is five dimensional. However, producing tables of as a function of five parameters is impractical. To reduce the task, we make the assumption that the corrected time PDF depends only on the track energy, the predicted “prompt” charge, and the predicted “late” charge, where “prompt” corresponds roughly to light arriving at the PMT directly from the track without delays induced by emission lifetimes, scattering, etc. Loosely, the assumption is that the shape of the corrected time spectrum is dominated by the physical extent of the track (characterized by its energy), and by the amount of prompt and late light reaching the PMT. The extent of the track affects the spread of possible hit times, since there is a spread of photon production times, while the amounts of prompt and late light affect the proportion of peaked “prompt” response to the tail of “late” response.
This assumption reduces the configuration to three dimensions. We make one further simplification by assuming that for a given energy , can be modeled by separate primitive prompt and late distributions, indexed by and , respectively, from which we can construct the full PDF. As a result, the PDF is indexed not by the fundamentally twodimensional space of , but by and separately, with the full PDF calculated onthefly as described below. In practice, Cherenkov and scintillation primitive distributions are created and used as proxies for the desired prompt and late distributions, respectively.
The Cherenkov (prompt) primitive distributions are created by simulating particles throughout the detector with isotropically chosen directions and fixed energy . The particles are created with direct Cherenkov light only; all other sources of light (scintillation, scattering, etc.) are turned off. For each event, the true track parameters are used to evaluate the direct Cherenkov predicted charge for each hit. Histograms of the corrected time of the hits are produced for various ranges of predicted charge. Figure 10 shows three such distributions for muons. Since only direct Cherenkov light is present, the histograms show no latetime features. The shapes of the time spectra depend on the values involved, with the distributions becoming narrower and earlier with increasing predicted charge due to the increasing probability that an early photon will be recorded at the PMT. The distribution in each predicted charge range is fit to a Gaussian distribution, and the resulting Gaussian parameters (mean and width) are subsequently parametrized across using a sixthorder polynomial in . Figure 11 shows an example of these “secondlevel” fits for muons. The procedure is repeated at many values of , with the two secondlevel fits providing seven parameters for the Gaussian mean and seven parameters for the Gaussian width at each energy. The energy dependences of these 14 parameters are then fit as a function of to fourthorder polynomials in a thirdlevel parametrization, as shown in Figure 12 for the first seven of the secondlevel parameters. In addition to conveniently summarizing the dependence of the time response on and , the parametrizations provide a smooth likelihood surface, as required by the minimization algorithm.
This completes the prompt primitive distributions . They are calculated for a given predicted charge at a PMT for a track with energy by evaluating the 14 thirdlevel curves which parametrize the dependence of the secondlevel functions. These 14 values then determine the two functions which parametrize the dependence of the mean and width of the Gaussian parameters on at . These are evaluated at the particular value of to determine the appropriate mean and width of the Gaussian distribution which describes the Cherenkov primitive distribution appropriate for the particular values of and . The validity of the assumptions and simplifications can be tested directly by comparing the calculated primitive distributions with those actually observed in the Monte Carlo simulated events. Figure 13 shows such a comparison.
For the scintillation (late) primitive distributions, events are generated with scintillation light only, and the histograms are created in bins of . Each histogram is fit to a sum of two exponentials convolved with a Gaussian distribution. The exponential decay constants are fixed at and , leaving three free parameters: the time origin, the Gaussian resolution, and the relative weight of the two exponentials. At a given , the three parameters are extracted in each histogram to obtain their dependence on , which is parametrized as a sixthorder polynomial in . The dependence of the three sets of seven parameters from the polynomials are parametrized as a function of by fourthorder polynomials, in analogy with the Cherenkov case. We now have , the scintillation primitive distribution. Figure 14 shows the primitive distributions compared with actual distributions from the Monte Carlo simulation for muons at , 600 and .
To obtain the complete PDF for a given PMT hit, we first divide the predicted charge into “prompt” and “late” components based on the sources:
(26) 
These definitions of and encapsulate the following assumptions. First, all prompt light is due to direct Cherenkov light; all other light (including indirect Cherenkov light) follows the late scintillationbased time distribution. The of direct Cherenkov light which is included in accounts for the late pulsing of PMTs which causes promptly arriving light to appear at late times. The late time distribution is used to model indirect Cherenkov light since the fluorescence, scattering, and reflection processes which lead to indirect light have a time structure similar to that of scintillation. Further, late timing is less critical to the reconstruction, as prompt light provides essentially all the useful timing information.
The and values are used to combine the prompt and late primitive distributions to form the total PDF . In the process, one must account for the fact that a PMT can register only one hit for a given track due to the latency period of the electronics. Since the time reported by the PMT reflects the time of the first photoelectron, we assume that a PMT hit that includes a prompt photoelectron obeys the prompt primitive distribution regardless of how many late photoelectrons may follow. The Poisson distribution gives us the probabilities that no prompt or late photoelectrons are produced based on their expected mean values:
(27) 
With these probabilities, one can calculate the probability that a given hit contains at least one prompt photoelectron:
(28) 
We assign this as the weight for the prompt primitive distribution within the total PDF, while is used as the weight for the late primitive distribution, yielding:
(29) 
where and are the prompt Cherenkovbased and late scintillationbased timing distributions, respectively. Note that Eq. (28) removes the overall probability of a hit occurring. Thus, even if the absolute probability of a prompt photoelectron is small, if .
5 The Twotrack Reconstruction
The single track model with seven parameters is sufficient for reconstructing and CCQE events as well as cosmic muons or their decay electrons. As discussed in Section 2, NC events require a twotrack model. Figure 15 shows the 12 parameters needed to describe two photon tracks originating from a common vertex. The electron track model serves doubleduty as the photon track model, under the assumption that a showering electron is indistinguishable from a showering photon apart from the conversion distance. (The mean conversion length in Marcol 7 is 67 cm.) While the twotrack model is conceptually a straightforward extension of the onetrack case, complexities of the likelihood space require special treatment in the likelihood maximization algorithm.
5.1 Twotrack Charge Likelihood
Since the charge likelihood depends only on the measured charge and the total predicted charge at each PMT, the twotrack charge likelihood can be obtained by adding together the predicted charges from the two tracks to obtain the total predicted charge at each PMT:
(30) 
5.2 Twotrack time likelihood
The onetrack time likelihoods account for the “first photoelectron only” nature of the electronics by calculating the probability that a prompt photoelectron is present and by then weighting the prompt primitive distribution’s influence on the full time PDF accordingly. This scheme is extended to handle two tracks as follows.
The two singletrack primitive distributions, and are formed for each track, where labels the track number. In anticipation of aggregating all late light later, the late “scintillation” primitive distributions are averaged to form:
(31) 
Of the two tracks, one will have a midpoint that is nearer to the target PMT than the other. This track’s quantities are labeled with “n” (near); the other’s, “f” (far). In analogy with Eq. (26), we define:
(32) 
The probabilities of not obtaining a photoelectron from these sources are:
(33) 
from which the following weights are obtained:
(34) 
where is the probability that a prompt photoelectron from the near track exists given that any photoelectron exists, etc. The weights are used to combine the two prompt primitive distributions and the averaged late primitive distribution:
(35) 
resulting in the complete twotrack corrected time PDF.
6 Minimization of
6.1 Onetrack
With the likelihood defined, the parameter set that minimizes its negative logarithm must be found. In the single track minimization process, two issues arise.
First, the energy parameter is tied to the geometry of the event via the track profiles . If the spatial parameters (, etc) are varied simultaneously with , the minimization algorithm can get confused by this correlation. The solution is to iterate the minimization process, as described below.
A second issue arises from the discrete PMT lattice, which imprints a small fluctuating signal on despite the smooth parameterizations of the input tables. Minimization algorithms that rely on gradients or that do not have controllable step sizes are thus poorly suited to the problem. Therefore, the SIMPLEX method in Minuit[19] is chosen over the MIGRAD method.
To initiate the minimization process, a vector of seed parameters is obtained from a fast fitter. In the first iteration of the minimization, the energy is held at its seeded value while the six remaining parameters (, , , , , ) are varied via Minuit/SIMPLEX to find a temporary minimum of . These six parameters are then fixed, while is freed for a second iteration. Finally, is fixed once more with the six other parameters varied (as in the initial iteration), to find the final minimum of F. The parameters from the final SIMPLEX call are returned as the bestfit set .
6.2 Twotrack
Figure 16 shows two situations that may arise in fitting a event. The event has two photon tracks whose Cherenkov rings are represented in black, while green and red represent two ring configurations corresponding to different starting values for the parameters. In case (a), the 12 parameters are near the correct configuration. In case (b), the parameters are such that both tracks are directed toward the dominant ring. The latter case corresponds to a local minimum which may offer a worse likelihood than case (a); however, the minimization algorithm may be trapped by the fact that any small changes to the parameters result in an increased value of . For example, sweeping one track over to the smaller black ring involves passing through a region with little detected Cherenkov light. The intermediate states form a barrier of disfavored values of that the minimization algorithm code will have difficulty overcoming, especially if several parameters must be adjusted simultaneously even to make the distant configuration an improvement.
The twotrack fits require a minimization approach that avoids these traps in the likelihood surface. Monte Carlo events were used to identify both common and rare trapping scenarios, and the minimization algorithm addresses each of these through the following collection of seed configurations and minimization sequences:n

The conversion lengths and are seeded with either 50 cm or 250 cm, leading to four possible pairs.

and are seeded with the results from the electron singletrack fit or with one of the eight perturbations (see below).

and are seeded with the “best” directions from a full grid of tested directions (see below).

The fourvertex of the event is seeded with the fourvertex from the electron singletrack fit shifted according to , and .

is seeded with approximately from the electron singletrack fit, while the seed for is, when possible, the energy needed to give an invariant mass equal to .
For the last item, the energy seeds are based on an empirical expression that accounts for the second ring’s energy contribution to the onetrack fit energy as a function of the angle between the two tracks. The function is obtained by simulating NC events and reconstructing them with the single track electron fit. Using the true photon energies () and opening angle , the relationship between and is parametrized by a linear function. This function is used to determine the seed and values for a given and .
The nine () possible seed directions for track 1 are created as follows. The spatial distribution of PMT charge is projected onto a plane perpendicular to the electron fit direction. The major axis of the covariance ellipse of the resulting twodimensional distribution is found. The nine possible seed directions are rotations of the best electronfit direction by , , and radians parallel to the major axis and radians perpendicular to the major axis.
The seed directions for track 2 come from either a “coarse” grid or a “fine” grid. The grid type is specified at run time.
Of the or possible combinations of track 1 and track 2 seed directions, only six are used to seed actual minimization sequences. They are the six combinations with (1) the best total likelihood (time and charge), (2) the best charge likelihood when the estimated track energies
The six best and secondbest track 1 and track 2 direction and energy combinations are chosen anew for each permutation of , leaving 24 complete parameter sets in total. Each set seeds two Minuit sequences. In sequence 1, the SIMPLEX minimization routine is first run with the energy and angle parameters held fixed; once a minimum is found, SIMPLEX is called again with all parameters free. In sequence 2, only the energy parameters are held fixed in the initial SIMPLEX call. This results in 48 fit sequences, each of which ends with a completely free SIMPLEX minimization. The best parameter set (i.e., the one that results in the lowest ) from the fortyeight final SIMPLEX calls is reported as the best parameter set for the twotrack fit. During all minimizations, the photon conversion points are constrained to be within the main detector volume. This prevents the fit from removing the influence of a track through an inflated conversion length. The event vertex itself, however, is not constrained.
The SIMPLEX minimizations can be performed with a constraint on the invariant mass. This is accomplished by removing as a free parameter and setting it via the relation:
(36) 
where is the angle between the photon tracks. This fixedmass mode is the actual hypothesis, while the freemass mode allows for mass reconstruction. Both are used for identification, with the former lending its maximized likelihood and the latter providing a reconstructed mass .
7 Performance
The performance of the reconstruction algorithm on simulated neutrino interactions
( CCQE, CCQE and NC events) generated with NUANCE [16] and
propagated through the detector Monte Carlo simulation is shown in Figures 1721.
A few observations about the figures, particularly Figure 18: (1) The resolutions for CCQE events degrade as the energy approaches the Cherenkov threshold from above, as most of the information in the event comes from the (rapidly decreasing) prompt Cherenkov light. (2) A lowenergy offers little information on its direction, as all that is observed is the isotropically produced final state. Thus, the angular resolution (but not the vertex resolution) becomes poor. (3) The vertex resolution for events is worse than CCQE events since the former has spatially displaced light production thanks to the average photon conversion length. (4) Typical performance numbers for CCQE events: 20 cm vertex resolution, kinetic energy resolution, and angular resolution. (5) Typical performance numbers for CCQE events: 10 cm vertex resolution, kinetic energy resolution and angular resolution. (6) Figure 19 shows a mass resolution of .
8 Particle Identification
The maximum likelihoods returned from each fit can be used for hypothesis testing. The two quantities
(37) 
and
(38) 
are used to determine for a given event whether the electron hypothesis is preferred over the muon and hypotheses. In these expressions, , and are the maximized likelihoods returned by the electron, muon, and (fixedmass) twotrack fits, respectively.
8.1 Electron/Muon separation
Figure 22 shows for simulated CCQE events and CCQE events the distribution of as a function of the energy reconstructed by the electron fit. For each sample, the events have been subject to a set of preselection criteria used in the appearance analysis. These require that there is only one timecluster of PMT hits in the event, eliminating obvious CC events that produce decay electrons. Cosmic backgrounds are eliminated by requiring the minimum number of PMT hits in the main region to be and the number of veto hits to be . Furthermore, the average time of hits is required to lie within the expected beam delivery window from the Booster.
For the CCQE events, tends to take positive values, indicating that the fit to the event with the electron hypothesis is favored over the muon hypothesis. Likewise, tends to be negative for the CCQE events, indicating that the muon hypothesis fits these events better than the electron hypothesis. At high energies (), the electron/muon separation is aided by the fact that the muon pathlength grows approximately linearly with energy, while the electron shower grows more slowly.
8.2 Electron/ separation
The left plots in Figure 23 show the distributions of as a function of energy reconstructed by the singletrack electron reconstruction for simulated CCQE events and NC events. The quantity tests whether a given event fits better as a single electron track or as a . Overall, electron/ separation becomes more difficult at high energies as the energies of the decay photons in events become more asymmetric and the shower fluctuations in a single electron event get large enough to mimic the presence of a second photon. Also, events in which one of the two photons leaves the detector unconverted present an essentially irreducible background at low energy.
Another quantity that can be used for electron/ separation is , the invariant mass of the two photons returned by the freemass twotrack fit. The NC events peak at the mass, whereas the CCQE events peak toward zero. This is seen in the right plots in Figure 23, where the distribution is shown as a function of reconstructed energy. The use of the true mass in the fit’s seeding procedure (Section 6.2) induces no appreciable highmass peak in the CCQE events while helping considerably in the correct identification of NC events.
The appearance analysis at MiniBooNE utilizes both and to separate CCQE events from NC .
8.3 Particle Identification Performance
The background rejection levels reachable with these fitterderived quantities are shown for Monte Carlo events in Figures 24 and 25. Figure 24 shows the efficiency, after fiducial volume and cosmic rejection cuts, for selecting signal CCQE events alongside the efficiency (or misidentification rate) for CCQE events. Three prototype selections are used, each chosen to give an approximately fixed signal efficiency at all energies. Note that the misidentification rate is shown with a logarithmic vertical scale. Figure 25 shows similar plots for rejection, one for the cut and one for the cut.
Misidentification rates of a few percent for CCQE events are seen for signal efficiencies from . (Note that CCQE events are further reduced by about an order of magnitude by recognizing the delayed muon decay.) The two variables indicate that the misidentification rate increases with energy; this is expected, since the two decay photons become highly asymmetric energetically and since the Cherenkov rings can achieve a greater degree of overlap. The rise in the misidentification rate at low energies is due to events with an escaping photon (corresponding to the cluster of events with positive at 50200 in Figure 23). Between 200 and 600 , where most of the NC events lie, a misidentification rate of less than 1 percent can be achieved with efficiency using only the cut, while rates of a few percent can be achieved with only the cuts. In the MiniBooNE oscillation analysis, 50% signal efficiency and 1% misidentification rate are obtained by requiring both and , where the “max” functions are quadratic in the energy obtained from the singletrack electron fit. These quadratic selection boundaries are indicated by the solid black lines in Figure 23. These cuts, along with an analogous cut, were tuned to optimize sensitivity to LSNDlike oscillations.
9 Summary
A maximum likelihood reconstruction algorithm is used at MiniBooNE to determine track parameters under electron, muon, and hypotheses. Underlying the likelihood is a track and detector model that calculates for a given parameter set the charge and time PDFs expected for each PMT, accounting for the spatially extended production of Cherenkov and scintillation light as well as the effects of indirect light from subsequent optical processes.
The maximized likelihoods for a given event obtained under different event hypotheses are used for event selection. In particular, the ratio of the likelihoods under the electron and muon models is used to suppress CC events in the MiniBooNE oscillation search; the ratio of the likelihoods under the electron and twotrack fit (with fixed invariant mass) and the invariant mass obtained from the freemass twotrack fit are used to suppress NC events. While the reconstruction has been developed within the context of the MiniBooNE oscillation search, other experiments employing similar Cherenkov detection techniques should be able to use the techniques discussed here.
Footnotes
 thanks: Now at the California Institute of Technology, Physics Department 10333, Pasadena, CA 91125, USA.
 thanks: Now at the University of British Columbia, Department of Physics and Astronomy, Vancouver, BC V6T 1Z1, Canada.
 Recall the energy seed procedure discussed above.
 For comparisons of Monte Carlo to data and for discussion of the development and validation of the Monte Carlo, see Refs. [17] and [18].
References
 kamioka
 K. S. Hirata et al. [KamiokandeII Collaboration], Phys. Rev. Lett. 65, 1301 (1990).
 Y. Fukuda et al. [Kamiokande Collaboration], Phys. Lett. B 335, 237 (1994).
 superk
 S. Fukuda et al. [SuperKamiokande Collaboration], Phys. Lett. B 539, 179 (2002) [arXiv:hepex/0205075].
 M. H. Ahn et al. [K2K Collaboration], Phys. Rev. Lett. 90, 041801 (2003) [arXiv:hepex/0212007].
 S. Yamamoto et al. [K2K Collaboration], Phys. Rev. Lett. 96, 181801 (2006) [arXiv:hepex/0603004].
 sno
 Q. R. Ahmad et al. [SNO Collaboration], Phys. Rev. Lett. 87, 071301 (2001) [arXiv:nuclex/0106015].
 Q. R. Ahmad et al. [SNO Collaboration], Phys. Rev. Lett. 89, 011301 (2002) [arXiv:nuclex/0204008].
 S. N. Ahmed et al. [SNO Collaboration], Phys. Rev. Lett. 92, 181301 (2004) [arXiv:nuclex/0309004].
 A. A. AguilarArevalo et al. [The MiniBooNE Collaboration], Phys. Rev. Lett. 98, 231801 (2007) [arXiv:0704.1500 [hepex]].
 A. A. AguilarArevalo et al. [The MiniBooNE Collaboration], Nucl. Instrum. Meth. A 599, 28 (2009).
 S. J. Brice et al., Nucl. Instrum. Meth. A 562, 97 (2006).
 CERN Program Library Long Writeup W5013 (1993).
 D. Casper, Nucl. Phys. Proc. Suppl. 112, 161 (2002) [arXiv:hepph/0208030].
 R. B. Patterson, Ph.D. thesis, Princeton University, FERMILABTHESIS200719 (2007)
 A. A. AguilarArevalo et al. [MiniBooNE Collaboration], Phys. Lett. B 664, 41 (2008) [arXiv:0803.3423 [hepex]].
 CERN Program Library Long Writeup D506