Measurement of the forward-backward asymmetry in top quark-antiquark production in ppbar collisions using the lepton+jets channel

We present a measurement of the forward--backward asymmetry in top quark-antiquark production using the full Tevatron Run II dataset collected by the D0 experiment at Fermilab. The measurement is performed in lepton+jets final states using a new kinematic fitting algorithm for events with four or more jets and a new partial reconstruction algorithm for events with only three jets. Corrected for detector acceptance and resolution effects, the asymmetry is evaluated to be 10.6+-3.0 %. Results are consistent with the standard model predictions which range from 5.0% to 8.8%. We also present the dependence of the asymmetry on the invariant mass of the top quark--antiquark system and the difference in rapidities of top quark and antiquark.


A. Motivation and definitions
Over the last five years both experiments at the Fermilab Tevatron Collider measured positive forwardbackward asymmetries in the production of top quarkantiquark pairs in proton-antiproton collisions (pp → tt) [1][2][3][4][5]. The reported values were consistently above predictions of the standard model of particle physics (SM) [6,7]. In particular, the CDF Collaboration observed a strong rise of the asymmetry with the invariant mass of the tt system, m tt [3]. The dependence of the asymmetry on m tt in D0 data, as measured in Ref. [4], was statistically compatible with both the SM predictions and with the CDF result. Several beyond-the-SM scenarios were suggested to explain the measured A FB values [8], in particular using the framework of parityviolating strong interactions suggested in Ref. [9]. In this paper we report new results from the D0 experiment based on the full dataset collected during Run II of the Fermilab Tevatron Collider, which supersede the result of Ref. [4].
In proton-antiproton collisions, top quark-antiquark pairs are predominantly produced via valence quarkantiquark annihilation. Thus, the direction of the proton (antiproton) almost always coincides with the direction of the incoming quark (antiquark). We define the difference in rapidity 1 between the top quark (y t ) and antiquark (yt): ∆y = y t − yt. (1) We refer to the events that have ∆y > 0 as "forward", and to those with ∆y < 0 as "backward". The forwardbackward asymmetry in tt production is defined as where N f (N b ) is the number of forward (backward) events. All tt asymmetries reported in this paper are given after subtracting the contributions from background processes. The rapidities of the t andt quarks and the corresponding asymmetries can be defined at the production level (sometimes denoted as generator level or parton level), when the kinematic parameters of the generated top quarks are used. Unless stated otherwise the productionlevel asymmetries are defined for all signal events without imposing the selection criteria of this analysis. The rapidities and asymmetry can also be defined at the reconstruction level, using the reconstructed kinematics of the selected events. Similarly, the invariant mass of the tt system can be defined at the production and reconstruction levels.

B. Strategy
In the SM a top quark almost always decays to a b quark and a W boson, which decays either leptonically or hadronically. In this paper we identify tt events using the tt → W + bW −b ; W + → l + ν l ; W − → qq ′ (and charge conjugates) decay chain. This channel is commonly referred to as the "lepton+jets" (l+jets) channel. We select events that contain one isolated lepton (electron or muon) of high transverse momentum (p T ) and at least three jets. The electric charge of the lepton identifies the electric charge of the leptonically decaying W boson and its parent top quark. The other top quark is assumed to have the opposite charge. The event selection, sample composition determination, and modeling of the signal and background processes are identical to those used in the measurement of the leptonic asymmetry in tt production in the l+jets channel [10]. The four-vectors of the top quarks and antiquarks in the events containing at least four jets are reconstructed with a kinematic fitting algorithm, while for the events that contain only three jets a partial reconstruction algorithm is used. If a jet exhibits properties consistent with a jet originating from a b quark, such as the presence of a reconstructed secondary vertex, we identify it as a b-tagged jet [11]. The l+jets events are separated into channels defined by jet and b-tag multiplicities. The amount of signal and the forward-backward asymmetry at the reconstruction level are determined using a simultaneous fit to a kinematic discriminant in these channels.
The measured background-subtracted one-dimensional (1D) distribution in ∆y is corrected to the production level ("unfolded"). From this distribution we calculate the fully-corrected A FB as well as its dependence on ∆y.
To study the dependence of the asymmetry on the invariant mass of the tt system, unfolding is done on the background-subtracted data distributions in two dimensions (2D: ∆y vs m tt ). The signal channels are unfolded simultaneously to yield the desired 1D or 2D productionlevel distributions, from which the production-level A FB values are computed using Eq. 2. The procedure is calibrated using simulated samples with varied asymmetries and input distributions in ∆y and m tt . The statistical and systematic uncertainties of the results are evaluated using ensembles of simulated pseudo-datasets (PDs).

II. D0 DETECTOR
We use the data collected by the D0 detector during Run II of the Tevatron in the years 2001-2011. After imposing event quality requirements ensuring that all detector systems were fully operational, this dataset corresponds to an integrated luminosity of 9.7 fb −1 . The D0 detector is described in detail elsewhere [12]. The central tracking system, consisting of a silicon microstrip tracker and a scintillating fiber tracker, is enclosed within a 1.9 T superconducting solenoid magnet. Tracks of charged particles are reconstructed within a detector pseudorapidity region 2 of |η| < 2.5. Electrons, photons, and jets of hadrons are identified [13] using a liquid-argon and uranium-plate calorimeter, which consists of a central barrel covering up to |η| ≈ 1.1, and two endcap sections that extend coverage to |η| ≈ 4.2 [14]. Central and forward preshower detectors are positioned in front of the corresponding sections of the calorimeter. A muon system consisting of layers of tracking detectors and scintillation counters placed in front of and behind 1.8 T iron toroids [15] identifies muons [16] within |η| < 2. Luminosity is measured using arrays of plastic scintillators located in front of the endcap calorimeter cryostats. A three-level trigger system selects interesting events at the rate of ≈ 200 Hz for offline analysis [17].

III. EVENT SELECTION AND MODELING
Object reconstruction and identification, as well as event selection, are the same as in Ref. [10] and are briefly outlined in this section. We select events with exactly one isolated electron within the detector pseudorapidity range of |η| < 1.1 or one isolated muon within |η| < 2.0, and at least three jets within |η| < 2.5. To limit the possible contribution of poorly modeled background due to multijet production, leptons of either flavor are required to have |y| < 1.5. The presence of a neutrino is inferred from a transverse momentum imbalance, which is measured primarily using calorimetry and is referred to as the "missing transverse energy" / E T . All selected objects are required to have transverse momentum p T > 20 GeV, and the jet with the largest p T (the leading jet) is also required to have p T > 40 GeV.
To identify jets that are likely to be associated with b quarks, we perform a multivariate analysis (MVA) that combines variables characterizing the properties of secondary vertices and of tracks with large impact parameters relative to the primary pp interaction vertex (PV) [11]. The output of the MVA is a continuous variable, MVA b . The requirement on MVA b (b tagging) used in this analysis has an efficiency of about 64% for identifying b jets originating from top quark decay, and a misidentification probability of about 7% for jets that do not contain heavy flavor quarks and are produced in association with leptonically decaying W bosons.
We simulate tt production using MC@NLO program (version 3.4) [18] with the parton showering performed by HERWIG [19]. This simulation is fully integrated with the D0 software, allowing for detailed studies of the kinematic dependences of A FB and their interplay with selection and reconstruction effects. The main source of background to the tt signal is the production of a leptonically decaying W boson in association with jets (W+jets). The kinematic properties of this process are simulated using ALPGEN [20] with hadronic showering performed by PYTHIA [21]. For signal and background modeling we use the CTEQ6.1 set of parton distribution functions (PDFs) [22]. The normalization of the W+jets contribution is a free parameter in the fitting procedure described below. Events with multiple jets can also mimic tt signal when a particle from one of the jets is misidentified as an isolated lepton. The normalization of this multijet background is extrapolated from a control sample enriched in this process using the probability for a jet to satisfy the lepton-quality requirements [23]. For the other backgrounds, Z+jets events are simulated with ALPGEN, diboson events are simulated with PYTHIA, and events from single-top-quark production are simulated with COMPHEP [24]. The normalizations for the last three background processes are taken from NLO calculations [25]. In all cases, event generation is followed by the GEANT-based D0 detector simulation [26]. To model energy depositions from noise and additional pp collisions within the same bunch crossing, simulated events are overlaid with data from random pp crossings. All simulated events are reconstructed using the same code as for the reconstruction of the collider data.

IV. RECONSTRUCTION OF THE EVENT KINEMATICS
To measure the tt forward-backward asymmetry and its dependence on the invariant mass of the tt system we need to determine the four-vectors of top quark and antiquark, which is done by summing the four-vectors of their decay products in the tt → W + bW −b ; W + → l + ν l ; W − → qq ′ (and charge conjugates) decay chain. There are four final state quarks in this decay chain, while we select events that contain one isolated lepton and at least three jets. When an event contains at least four jets we assume that the four jets with the largest p T originate from the quarks from tt decay. If an event contains only three jets one of the jets from tt decay is missing. In either case all possible assignments of three or four jets to the final state quarks are used, with the likelihood of each assignment evaluated by the reconstruction algorithms described below.
For l+≥4 jet events, the tt system is fully reconstructed using a kinematic fitting algorithm. Previous D0 top quark analyses used the algorithm of Ref. [27]. In this paper a new algorithm is employed, which utilizes an analytic solution for the neutrino momentum using the constraints on the W -boson (M W ) and top-quark masses (m t ) [28]. The likelihood term for each jet-to-quark assignment accounts for the differences between the observed jet energy and the energy scaled to satisfy the constraints on M W and m t . The jet energy resolution and the probability for a jet to be reconstructed (see "Type III" transfer function in Ref. [29]) are taken into account. The b-tagging observables MVA b are also used to evaluate the likelihood of each assignment.
For l+3 jet events, a partial reconstruction algorithm of the tt decay chain is employed [30]. With one jet entirely lost, no significant improvement is expected from scaling the four-vectors of the remaining objects as is done by the kinematic fitting algorithm in l+≥4 jet events, so the partial reconstruction algorithm does not attempt to modify the kinematics of the observed objects. As only the transverse components of the neutrino momentum are measured in / E T , the longitudinal component is calculated using a quadratic equation which results from imposing the M W constraint on the W → lν decay products. The two-fold ambiguity is resolved by choosing the solution that minimizes the difference between the known m t and the invariant mass of the objects assigned to the leptonic top quark decay, m l . This algorithm thus assumes that the jet associated with the b quark from the leptonically decaying top quark is detected. This assumption holds for 80% of the tt events. The lost jet is assumed to be associated with either a light quark or a b quark from the hadronic top-quark decay chain. In the majority of cases (74%) this jet is lost due to its low energy, so this loss has little effect on the kinematics of the hadronically decaying top quark. The lost jet is neglected in the partial reconstruction algorithm. The sum of the four-vectors of the two jets assigned to the products of the hadronically decaying top quark serves as a proxy for the four-vector of the hadronically decaying top quark with the invariant mass m p . Even though m p is not expected to be equal to m t , the distribution in this variable is different for combinations correctly associated with the hadronically decaying top quark and combinations that include a b jet from the leptonically decaying top quark. In each event we consider the following nine observables: the MVA b for each of the three jets, the three possible m l , corresponding to the three possible lepton-neutrino-jet combinations, and the three possible m p . The likelihood of each of the three possible jet-to-quark assignment is calculated by evaluating the consistency of the nine observables with the distributions corresponding to the hypothesized assignment. In particular, the jet hypothesized to be associated with a b quark should have a value of MVA b consistent with the one expected for b jets, while for a jet hypothesized to originate from a W boson decay MVA b should be consistent with the distribution expected for such jets. The values of m l and m p for the jet combinations that correspond to the hypothesized assignment should be consistent with the distributions expected for correctly assigned jets, while the values of m l and m p for the other jet combinations should agree with the distributions expected for wrong assignments. When calculating m tt , we compensate for the effect of the lost jet by applying an m p -dependent scaling to the four-vector of the hadronically decaying top quark.
Unlike the A FB measurement in Ref. [4], where only the jet-to-quark assignment with the lowest χ 2 was used, we reconstruct ∆y by averaging its values over all possible assignments, weighted by their likelihoods evaluated as described above for l+≥4 jet and l+3 jet events. The same approach is used to reconstruct m tt in the l+3 jet channel. For l+≥4 jet events, m tt is reconstructed using the outputs of three reconstruction algorithms: the new kinematic fit algorithm, the kinematic fit algorithm of Ref. [27], and a simple reconstruction algorithm [31] that evaluates the kinematics of the leptonically decaying W boson from the lepton and the neutrino by imposing the M W constraint and calculates m tt by adding the four most energetic jets without imposing the m t constraint. The likelihood values calculated by the algorithms give indications on how well the kinematics of a particular event match the assumptions made by a given algorithm. In particular, for high m tt there is a higher probability that two final state quarks are associated with the same jet. Such a jet is likely to be the most energetic jet in the event and have a large mass. The simple reconstruction algorithm, which does not assume a specific jet-to-quark assignment, performs best for such events. We use a multivariate regression [32] to combine the partially correlated m tt values and the likelihoods produced by the three algorithms with supplementary observables such as the mass of the leading jet to estimate m tt . This combined m tt reconstruction outperforms the individual algorithms in all m tt ranges.
For the asymmetry measurement the performance of a tt reconstruction algorithm can be characterized by the probability P c to correctly reconstruct the sign of ∆y. For the algorithm employed in this analysis for l+≥4 jet events P c = 0.775, compared to P c = 0.756 for the algorithm of Ref. [27]. The partial reconstruction algorithm achieves P c = 0.745 for l+3 jet events. The dependence of P c on the production-level |∆y| is shown in Fig. 1 for these three algorithms. The high values of P c achieved by the partial reconstruction algorithm, which are almost as high as P c for l+≥4 jet events, can be understood from the following consideration. All four leading jets are associated with the quarks from the tt decay in only 55% of the l+≥4 jet events. For the other 45% of the events one of the jets originates from initial or final state radiation, which can lead to badly misreconstructed tt four-vectors. Only 4% of the l+3 jet events contain a jet that does not originate from the four quarks of the tt decay. Thus, even though some information is lost with the unreconstructed jet, no wrong information is added, leading to a low probability to misreconstruct the sign of ∆y. . The probability to correctly reconstruct the sign of ∆y as a function of the production-level |∆y| for the algorithm of Ref. [27] used to measure the AFB in Ref. [4] and the algorithms used to reconstruct l+≥4 jet events and to partially reconstruct l+3 jet events in this paper.

V. SM PREDICTIONS FOR AFB
The differential cross section for tt production is available only at order α 3 s , where α s is the strong coupling constant. Since in the SM the tt asymmetry only appears at this order, no full higher order prediction for A FB exists yet. The relative uncertainty on the α 3 s calculation of the asymmetry due to higher order corrections is evaluated to be as large as ≈ 25% [33].
Recently the order α 4 s calculation for the total cross section of tt production [34] was made available, but the asymmetry was not computed at this order. Several papers report calculations of the leading corrections to the asymmetry with the predicted A FB values ranging from 5.0% in MC@NLO to 8.8% once the electroweak corrections [35] and resummations of particular regions of phase space [36] are taken into account. The dominant uncertainty on these predictions is from the renormalization and factorization scales, and is evaluated to be as high as 2.0% (absolute) [33,37]. The authors of Ref. [38] obtain a value of A FB = 12.7% by choosing a normalization scale that arguably stabilizes the perturbative expansion yet differs significantly from the scales commonly used in top quark physics calculations. Some authors suggest that the corrections from interactions between the top quark decay products and the proton remnants should also be taken into account when calculating A FB [39]. Given this variety of predictions, we choose to compare our data to the well defined MC@NLO simulation.
At order α 3 s , the QCD contributions to the asymmetry in tt production can be divided into two classes up to divergences that cancel between these two classes [6]. The first class, which contributes to negative asymmetry, is a result of interference between the terms that contain gluon radiation in the initial or final states, which may result in an extra jet in the event and typically leads to a higher transverse momentum of the tt system. The second class, which contributes to positive asymmetry, is from interference between the Born term (α 2 s ) and the term described by a box diagram (α 4 s ). The overall asymmetry is positive and depends on the jet multiplicity. Selection criteria that give preference to events with higher jet multiplicity favor the first class of events and further lower the overall expected asymmetry, while a higher asymmetry is expected for events with lower jet multiplicity. Consequently, forward events tend to have fewer jets than events in the backward category. Similarly, since a b-tagged jet is less likely to originate from initial or final state radiation, samples with a larger number of b tags tend to have higher values of A FB .
MC@NLO predicts an overall asymmetry in tt production before selection of (5.01 ± 0.03)%. Here, and in the following sections, the quoted uncertainties on the predictions are from the finite size of the simulated samples unless otherwise stated. Table I lists the MC@NLO predictions for tt events after the selection criteria are applied.
All previous measurements of A FB in the l+jets channel selected tt events that had at least four jets in the final state. As is apparent from Table I, restricting the selection to only l+≥4 jet events lowers the productionlevel asymmetry. Including events with three jets reduces this selection bias.
Asymmetries after reconstruction are presented in the last column of Table I. Finite resolution in ∆y results in roughly 20% of the forward events being misreconstructed as backward, and vice versa. Since there are more forward events, ∆y smearing leads to an overall lowering of the reconstructed asymmetries. At the same time, forward tt events, which tend to have fewer jets, have a lower probability to be misreconstructed, resulting in fewer migrations into the backward category, and an upward shift in the reconstructed asymmetry. This bias is most apparent in the l+≥4 jet, one-b-tag channel, where the lowest asymmetry is predicted.

VI. SAMPLE COMPOSITION AND RECONSTRUCTION-LEVEL AFB
Reconstructed events are divided into six channels by the number of jets and b tags: l+3 jet and l+≥4 jet with 0, 1, and ≥2 b tags each. The l+3 jet zero-b-tag channel is used only for the background asymmetry calibration, and not for the tt asymmetry measurement. The l+≥4 jet zero-b-tag channel is used only for determining the sample composition and the reconstruction-level A FB , and is not used for measuring the production-level asymmetry.
Several well-modeled variables that have different distributions for signal and background processes, and that have minimal correlations between each other and with ∆y and m tt , are combined into kinematic discriminants bounded between 0 and 1 [10]. For l+≥4 jet events a discriminant D 4 is built from the following input variables: • χ 2 -the test statistic of the likeliest assignment from the kinematic fit.
• p LB T -the transverse momentum of the leading btagged jet, or when no jets are b tagged, the p T of the leading jet.
between the two closest jets, a and b, and p T,a and p T,b are their transverse momenta.
• M jj , the invariant mass of the jets corresponding to the W → qq ′ decay in the likeliest assignment from the kinematic fit, calculated using kinematic quantities before the fit.
The variables χ 2 and M jj are based on the full tt reconstruction using the kinematic fitting technique of Ref. [27]. For the l+3 jet events we construct a discriminant D 3 using a different set of input variables: • S -the sphericity [40], defined as S = 3 2 (λ 2 + λ 3 ), where λ 2 and λ 3 are the two largest of the three eigenvalues of the normalized quadratic momentum tensor M . The tensor M is defined as where p o is the momentum vector of a reconstructed object o, and i and j run over the three indices for the Cartesian coordinates. The sum over objects includes the three selected jets and the selected charged lepton.  Fig. 2). The sum of the estimated number of signal and background events is constrained to be equal to that in data. The "Selected events" column includes the l+3 jet events with at least one b tag and all l+≥4 jet events. The statistical uncertainties from the fit are quoted. We also present the event breakdown for the channels with at least one b tag, which are used to determine the production-level AFB. • p 3rd T -the transverse momentum of the third leading jet.
• M min jj -the lowest of the invariant masses of two jets, out of the three possible jet pairings.
• p LB T , defined as for the l+≥4 jet channel, above. • ∆φ(jet 1 , / E T ), the difference in azimuthal angle between the leading jet and / E T .
The discriminants for all channels are combined into a single discriminant D c , so that for the l+3 jet events D c = N tag + D 3 , while for l+≥4 jet events D c = 3 + N tag + D 4 . The variable N tag above is usually taken to be equal to the number of b-tagged jets in the event, but for events with more than two b-tagged jets N tag = 2 instead. We fit the sum of the signal and background templates to the data distribution in the discriminant D c as shown in Fig. 2. This fit is identical to the fit for the sample composition in Ref. [10]. The sample composition and its breakdown into individual channels are summarized in Table II. Background contributions other than W+jets and multijet production are labeled "Other Bg" in Table II. In the simulated W+jets background, the angular distribution of leptons from W -boson decay has a forwardbackward asymmetry, which is in part tuned to Tevatron data [41]. Due to this asymmetry, when these events are reconstructed according to the tt hypothesis, there remains a residual asymmetry of ≈ 5% in the ∆y distribution. To improve the modeling of this asymmetry, we apply a weight to each simulated W+jets event which depends on the product of the generated lepton charge and its rapidity. These weights are chosen so that the simulation best matches control data with three jets and zero b tags as in Ref. [10]. The difference in the ∆y distributions predicted by the simulation with and without the applied weights is treated as a source of systematic uncertainty due to background modeling. This uncertainty exceeds the uncertainty due to PDFs by about a factor of two. We rely on the simulation to predict the variation of the asymmetry in W+jets events with jet and b-tag multiplicities.
The distributions of the reconstructed ∆y are shown in Fig. 3. The tt asymmetry at the reconstruction level is extracted using a fit to the distributions in the discriminant D c and sign of ∆y, excluding the l+3 jet events with zero b tags. This fitting procedure is identical to the procedure used in Ref. [10]. The inclusive asymmetry measured at the reconstruction level is (7.9 ± 2.3) %. The results for individual channels are listed in Table III.   TABLE III: Reconstruction-level background-subtracted asymmetries for selected events for different channels. The last row includes the channels listed in the rows above and the l+≥4 jet, zero-b-tag channel. The first uncertainty is statistical, and the second one is systematic. Systematic uncertainties are discussed in Section VIII. The prediction is based on the MC@NLO simulation. The distributions of the reconstructed invariant mass of the tt system are shown in Fig. 4. Since the l+3 jet and l+≥4 jet channels have different response (both mean and shape) for m tt , the dependence of A FB on m tt at the reconstruction level is difficult to interpret and is not presented here. The measurement of production-level A FB and its dependence on m tt is described in Section VII.
We use the results of the sample composition study summarized in Table II to normalize the distributions for the background processes in the sensitive variables (∆y, and also m tt for the 2D measurement), which are subtracted from the distributions observed in data. To increase the signal purity of the data used in the fully corrected measurements, we unfold only events containing at least one b tag. The background-subtracted distributions of the sensitive variables in the corresponding four channels are used as inputs to the unfolding procedure.

VII. UNFOLDING THE ASYMMETRY
The true or generated distribution of a certain variable (∆y for the inclusive measurement) is shaped by acceptance and detector resolution, resulting in the observed distribution, which is also subject to statistical fluctuations. The goal of the unfolding procedure is to find the best estimator for the true distribution given the background-subtracted data and knowing detector acceptance and resolution from simulation. After finding the  best estimator for the true distribution of ∆y, we summarize it into the production-level A FB using Eq. 2. This is the same general approach used in the previous measurement [4]. For this unfolding we use TUNFOLD [42], which we extend as discussed below.
Each distribution is presented as event counts in a binned histogram 4 , i.e., as a vector with a dimension equal to the number of bins. Given the vector of production-level tt signal counts, p, and the vector of expected background counts, b, the expected data counts in the i-th bind i is given bỹ where A is a diagonal acceptance matrix, whose jj-th element is the probability for an event produced in the j-th bin to pass the selection criteria and M is the normalized migration matrix, whose ij-th element is the probability for a selected event produced in the j-th bin to be observed in the i-th bin. Given the vector of observed counts, d, we can construct the vector of background-subtracted reconstruction-level counts, r = d − b, with its covariance error matrix V . The matrix V is constructed to account for the expected statistical uncertainties on data and background, in particular those due to the size of the multijetenriched control sample. We then seek to find the vector u, which best estimates the vector of production-level counts p, by minimizing for a given vector r, where τ is the regularization strength and R is the regularization matrix. The first term of Eq. 6 quantifies the consistency of u with data, while the second (regularization) term quantifies the smoothness of u.
Without regularization, the unfolding procedure amounts to a minimization of the first term in Eq. 6. If the numbers of reconstruction and production-level bins are equal, the problem of minimization is solved by simply inverting the matrix: u unregularized = T −1 r.
Unregularized matrix inversion typically results in unphysical, rapidly varying distributions [43]. Such distributions are disfavored in regularized unfolding by adding a second "regularization" term to the χ 2 . The regularization term in Eq. 6 depends on the discrete second derivative of the binned distribution u. For constant bin widths, the regularization term is calculated using a regularization matrix with the following structure [42]: For this analysis we modify the structure of R to regularize based on the second derivative of the event density 4 Overflows are included in the edge bins. rather than the event counts, which allows for the use of variable bin sizes. The regularization strength τ is chosen using both ensemble testing (described below) and the L-curve technique [42] to balance the minimization of statistical fluctuations and bias. The difference between the two techniques is included in the evaluation of the systematic uncertainty due to the choice of the regularization strength.
As in Ref. [4], the production-level ∆y distribution is divided into 26 bins and the reconstruction-level ∆y distribution is divided into 50 bins. Both have narrower bins near ∆y = 0, where the probability to misclassify forward events as backward or vice versa changes rapidly, and wider bins at high |∆y|, where statistics are low.
For the 2D measurement, we use six m tt bins at the production level, with edges at 0, 400, 450, 500, 550, 650 GeV and +∞. The joint distribution of ∆y and m tt has a kinematic boundary at |∆y| = log ( A bin edge close to this boundary would result in a large difference in the event density between adjacent bins, a feature that would be smoothed by a regularization procedure, thus biasing u. To avoid such a bias, the ∆y edges of the bins of the 2D measurement are chosen to depend on m tt as shown in Fig 5. The reconstruction-level histograms have similar but finer bins along both the ∆y and m tt directions. In the l+3 jet channels 13 m tt bins are used to accurately describe migrations among the six production-level bins. The m tt resolution in the l+≥4 jet channels allows for 14 m tt bins. We simultaneously unfold to the production level the four channels that contain at least one b tag. The difference in purity among channels is accounted for in the definition of the covariance error matrix V . The unfolding technique is calibrated, and the statistical and systematical uncertainties are determined using the results of ensemble tests. Each ensemble comprises simulated PDs that we build according to MC@NLO, ALPGEN [20] or MADGRAPH [44] SM predictions, or according to toy models with different asymmetries. The PDs are created from the expected bin countsd i calculated using Eq. 4 by adding Poisson (statistical) and Gaussian (systematic) fluctuations, with the Gaussian width taken as one standard deviation for the corresponding systematic uncertainty.
In the toy models the input distribution P (∆y) has the form: where a and δ are shaping parameters, G is a Gaussian distribution with mean µ and width wσ 0 , σ 0 the width predicted by MC@NLO, and w a scaling parameter. The shape of the ∆y distribution and the input asymmetry are varied using the parameters µ, a, δ, and w. In addition, we produce ensembles with the signal taken from simulated samples of tt production mediated by axigluons, hypothetical massive particles that arise in extensions of the SM that suggest different strong couplings for left and right-handed quarks [9]. The input asymmetry in the models used for calibration ranges from −30% to +30%, while the axigluon masses are varied from 0.2 to 2 TeV. The bias, which is the average difference between the unfolded and input A FB values, is shown in Fig. 6 as a function of the input A FB . Based on this study we derive a correction (calibration) that is applied to the result to eliminate the expected bias. The majority of the tested models are contained within the systematic uncertainty assigned to this calibration, which is shown by the dotted lines in Fig. 6. The one point that is significantly outside of the boundaries of this region corresponds to tt production mediated by an axigluon with a mass of 0.4 TeV. This particular model exhibits a significant change in A FB on the m tt scale smaller than the bin width (here, 50 GeV), thus breaking the assumption of a smooth underlying distribution, leading to biased results. The unfolding of models with such rapidly changing A FB will, in general, be biased in all regularized unfolding procedures, and we choose not to assign a systematic uncertainty that covers this specific class of models.
The unfolded ∆y distribution is presented in eight bins, with each bin calibrated for the expected bias observed in the MC@NLO-simulated PDs. The A FB value in each |∆y| range, A FB (|∆y|), is calibrated using the same procedure as for the inclusive A FB . Since no systematic correlation is found between the A FB biases in different |∆y|, as well as m tt bins, they are calibrated individually.

VIII. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties on the reconstruction and production-level A FB are summarized in Table IV in several categories, which are detailed below. To evaluate the systematic uncertainty on the reconstruction-level A FB , we vary the modeling according to the estimated uncertainty in the relevant parameter of the model and propagate the effect to the result. The systematic uncertainties on the production-level A FB are evaluated by including the effects of systematic variations on the simulated background-subtracted PD into the ensemble tests. To find the expected uncertainty due to each category we use dedicated ensembles generated without statistical fluctuations and with only the relevant systematic effects. The total uncertainties on the production-level A FB are taken from ensembles built including both statistical fluctuations and systematic effects (see Sec. VII). The background model category includes the following sources, which affect the properties predicted for background events. The leptonic asymmetry of the W+jets background is varied within its uncertainty of 3% [10]. The rate of heavy-flavor production within W+jets production is varied by ±20% [25,45]. The efficiencies for lepton identification, and the probabilities for a jet to be misidentified as a lepton, taken as functions of lepton momentum, are varied within their uncertainties to account for the uncertainty on the number of background events from multijet production [23]. This variation affects both the background shape and normalization. Uncertainties associated with the modeling of the discriminant, D c , transverse momentum of W boson and M min jj , as well as potentially increased background levels at high lepton pseudorapidity are also quantified by modifying the background model [10].
The signal model category includes the sources of uncertainty that affect the properties predicted for signal events other than the ones accounted for in the PDFs and pileup category. The top quark mass is varied according to the combined Tevatron measurement of Ref. [46]. The effect of higher order corrections to tt production is estimated by replacing the migration matrix M from Eq. 4 simulated by MC@NLO with the one simulated by ALPGEN, which uses tree-level matrix elements. The b quark fragmentation function is varied within its uncertainties [46], which also affects background modeling.
The signal model category also includes the uncertainties associated with gluon radiation. The total amount of initial state radiation is varied in a range consistent with the results of Ref. [47]. We also consider the difference in the predicted amount of initial state radiation between forward and backward events, both because of contributions at order α 3 s and due to higher order effects which are modeled by the simulated parton showers [48]. We account for this uncertainty by reducing the difference in the distributions of the p T of the tt system for forward and backward events by 25%, a value derived from Ref. [48]. We also account for the possibility that the mismodeling of this variable in the l+3 jet final state affects A FB by reweighting this distribution to match the D0 data, similarly to the procedure used in Ref. [10].
The uncertainties due to unfolding are dominated by the calibration uncertainties. The uncertainties associated with the choice of the regularization strength and statistical fluctuations in the MC samples used to find the migration matrix are also included.
The PDFs and pileup category includes uncertainties on the modeling of the pp collisions. The main uncertainties are from the PDFs, which primarily affect the ∆y distribution of the W+jets background. These uncertainties are evaluated by varying the contributions of the various eigenvectors from the CTEQ6.1 PDF [22] and by considering an alternative set of PDFs (MRST2003 [49]). The number of additional collisions within the same bunch crossing (pileup) affects the quality of the reconstruction. The uncertainties on the modeling of additional collisions are also included in this category.
The detector model category includes the following sources of systematic uncertainty. The efficiencies of the b-tagging algorithm for jets of different flavors, which are measured from collider data, are varied according to their uncertainties [11]. These variations affect the measured A FB mostly through the estimated sample composition, which depends strongly on the classification of data into several channels according to the number of b tags. The modeling of jet energy reconstruction, including the overall energy scale and the energy resolution, as well as jetreconstruction efficiencies and single-particle responses, are all calibrated to collider data and are varied according to their uncertainties [50]. The uncertainties due to jet reconstruction and energy measurement are significantly reduced compared to the previous measurement due to the inclusion of the l+3 jet events.
Lastly, the sample composition is varied according to its fitted uncertainties. This variation is performed in addition to the changes in the sample composition implicitly induced by other systematic variations.

A. Inclusive AFB and AFB dependence on |∆y|
The calibrated production-level ∆y distribution is shown in Fig. 7. The corresponding inclusive forwardbackward asymmetry in tt production is (10.6 ± 3.0) %. The dependence of A FB on |∆y| is shown in Fig. 8 and Table V with the correlation factors between bins listed in Table VI. These correlations are taken into account in the fit of the measured A FB (|∆y|) to a line. Since for any physical ∆y distribution the asymmetry at ∆y = 0 is 0, we constrain the line to the origin and fit for its slope. For data, we find a slope of 0.154 ± 0.043. This slope is compatible within two standard deviations with the MC@NLOsimulated slope of 0.080, which has negligible statistical uncertainty. The difference between the slope reported by the CDF Collaboration [5] and the slope reported in   The dependence of A FB on m tt is shown in Fig. 9 and Table VII with the correlation factors between bins listed  in Table VIII.  The values of the asymmetry measured in six m tt ranges constitute a six-dimensional vector v with a 6 × 6 covariance matrix Σ. Table IX lists the eigenvectors e i (i = 1, ..6) of Σ together with the corresponding components of the vector v in the basis formed by the eigenvectors: v i = v · e i , and their uncertainties σ i = Σ ′ ii , where Σ ′ is the covariance matrix transformed to the basis e i . The elements of Table IX fully specify the measured sixdimensional likelihood in the Gaussian approximation, and can be used for quantitative comparison with theoretical predictions and other experimental results [51].
Using the full covariance matrix we perform a fit of the measured A FB to the functional form We choose C = 445 so that the correlation factor between the fit parameters α and A 0 is less than 0.01 in the fit to the data. The parameters of the fit are listed in Table X for the data and the MC@NLO simulation. We observe a slope α consistent with zero and with the MC@NLO prediction. The difference between slope reported by the CDF Collaboration [5] and the slope reported in this paper corresponds to 1.8 standard deviations.  The dependence of the forwardbackward asymmetry on the invariant mass of the tt system. The D0 data points are shown with the total error bars indicating the total uncertainty, based on the diagonal elements of the covariance matrix, while the statistical uncertainties are indicated by the inner error bars. The dashed line shows a linear fit to the data with the dotted curves indicating the fit uncertainty. The horizontal lines correspond to the MC@NLO prediction for the asymmetry in each m tt bin [18]. Shaded boxes correspond to the prediction of Refs. [33,35]. The last bin has no upper boundary. The x coordinate of each data point is the simulated average of the m tt distribution in the corresponding bin.

X. DISCUSSION
The measured inclusive forward-backward asymmetry in tt production, A FB = (10.6 ± 3.0) % is in agreement with the SM predictions reviewed in Section V, which range from an inclusive asymmetry of 5.0% predicted by the MC@NLO simulation to (8.8 ± 0.9)% [35] once electroweak effects are taken into account. The measured dependences of the asymmetry on |∆y| and m tt are also in agreement with the SM predictions. Nevertheless, the observed A FB and the dependences of A FB on m tt and |∆y| do not disfavor the larger asymmetries that were previously measured in pp collisions [5].
To compare the presented result with the previous D0 publication [4], Table XI presents A FB at the reconstruction level measured in different samples. The method discussed in this paper applied to l+≥4 jet events from the first 5.4 fb −1 of integrated luminosity yields a result consistent with that in Ref. [4], but with a reduced uncertainty mainly due to the separation of data into channels based on the number of b tags and the increased efficiency of the new b-tagging algorithm. Once the analysis is extended to include the l+3 jet events collected at that time, the uncertainty is reduced by a factor of 1.26. The result obtained in the second 4.3 fb −1 of the Tevatron dataset is within one standard deviation from that obtained in the first 5.4 fb −1 . The statistical uncertainty obtained in the combined 9.7 fb −1 of integrated luminosity is reduced by a factor of 1.29 with respect to the result obtained using the same method in the first 5.4 fb −1 , while the reduction expected from scaling with the integrated luminosity is 1.34. This loss of sensitivity is mainly due to higher instantaneous luminosity during the collection of the later data, which required a tighter trigger selection. The improved reconstruction of ∆y and the reduced acceptance bias due to the inclusion of the l+3 jet events result in further reduction of the statistical uncertainty on the unfolded result compared to Ref. [4]. The separation of the data into channels allows us to add the l+3 jet channels without losing the statistical power of the purer l+≥4 jet channels.

XI. SUMMARY
In summary, we report the measurement of the forward-backward asymmetry in tt production using the dataset recorded by the D0 detector in Run II of the Fermilab Tevatron Collider. The results presented here supersede the ones that were based on about half of the data [4]. The analysis is extended to include events with three jets, allowing for the loss of one jet from the tt decay and reducing the acceptance corrections. The unfolding procedure now accounts for the differences in sample compositions between channels, thus maximizing the statistical strength of the individual channels. New reconstruction techniques are used in the l+≥4 jet channel, improving the experimental resolution in all variables of interest.
The asymmetry measured at the reconstruction level is A FB = (7.9 ± 2.3) %. After correcting for detector resolution and acceptance, we obtain a production-level asymmetry A FB = (10.6 ± 3.0) %. The observed asymmetry and the dependences of A FB on m tt and |∆y| are consistent with the standard model predictions.