Study of the K+ pi+ pi- Final State in B+ ->J/psi K+ pi+ pi- and B+ ->psi-prime K+ pi+ pi-

Using 535 x 10^6 B-meson pairs collected by the Belle detector at the KEKB e+ e- collider, we measure branching fractions of (7.16 +/- 0.10(stat) +/- 0.60(syst)) x 10^-4 for B+ ->J/psi K+ pi+ pi- and (4.31 +/- 0.20(stat) +/- 0.50(syst)) x 10^-4 for B+ ->psi-prime K+ pi+ pi-. We perform amplitude analyses to determine the resonant structure of the K+ pi+ pi- final state in B+ ->J/psi K+ pi+ pi- and B+ ->psi-prime K+ pi+ pi- and find that the K_1(1270) is a prominent component of both decay modes. There is significant interference among the different intermediate states, which leads, in particular, to a striking distortion of the rho line shape due to the omega. Based on the results of the fit to the B+ ->J/psi K+ pi+ pi- data, the relative decay fractions of the K_1(1270) to K rho, K omega, and K^*(892) pi are consistent with previous measurements, but the decay fraction to K_0^*(1430) is significantly smaller. Finally, by floating the mass and width of the K_1(1270) in an additional fit of the B+ ->J/psi K+ pi+ pi- data, we measure a mass of (1248.1 +/- 3.3(stat) +/- 1.4(syst)) MeV/c^2 and a width of (119.5 +/- 5.2(stat) +/- 6.7(syst)) MeV/c^2 for the K_1(1270).


I. INTRODUCTION
The large number of B-meson decays observed at B factories allows detailed studies of the intermediate-state resonances involved in these decays. This paper analyzes the structure of the K + π + π − final state in the decays B + → J/ψK + π + π − and B + → ψ ′ K + π + π − . 1 Kaon excitations that decay to a Kππ final state are difficult to distinguish based on the mass of the Kππ system alone, owing to their overlapping line shapes. 2 In this analysis, data are therefore fitted in the three dimensions M 2 (Kππ), M 2 (Kπ), and M 2 (ππ), which are the squared invariant masses of the K + π + π − , K + π − and π + π − systems, respectively. An unbinned maximumlikelihood fit is performed to extract maximal information from the data. The fitting model accounts for interferences among different intermediate states, as well as the spin-dependent angular distributions of the final state. The large sample size, combined with the clean environment afforded by the presence of a J/ψ or ψ ′ in the final state, makes it possible to distinguish the different kaon excitations that contribute to the K + π + π − final state. The results provide information not only on intermediate-state interactions but also on the structure of the kaon spectrum. By performing an additional fit in which the mass and width of the K 1 (1270) are floated, we measure the mass and width of the K 1 (1270). 1 Charge-conjugate modes are always implicit. 2 In 2001, the Belle Collaboration measured the branching fraction for B + → J/ψK 1 (1270) with 2% of the data presented here. The K + π + π − final state in B + → J/ψK + π + π − was found to be dominated by the K 1 (1270), and no other structure was detected [1].
Identifying the kaon excitations involved in B + → J/ψK + π + π − and B + → ψ ′ K + π + π − can lead to a better understanding of the underlying theory. For example, the breaking of SU (3) flavor symmetry mixes the 1 3 P 1 and 1 1 P 1 states of the kaon system into the physical states K 1 (1270) and K 1 (1400) as tromagnetic calorimeter, comprising 8736 CsI(Tl) crystals, records the energy deposited by photons, leptons, and hadrons. These subdetectors are surrounded by a superconducting solenoid 3.4 m in diameter and 4.4 m in length, which produces a 1.5-T magnetic field parallel to the positron beam. An iron flux return installed outside the coil is instrumented with large-area resistiveplate counters to identify muons and K L mesons. Monte Carlo (MC) simulations 3 are used to determine the acceptance of the detector for the processes of interest.

III. EVENT SELECTION
Electron candidates are identified by combining information from the CDC, electromagnetic calorimeter, and ACC. Muon candidates are identified by extrapolating charged-particle tracks from the silicon vertex detector and CDC into the K L /µ detector. To identify charged hadrons, momentum measurements from the CDC are combined with velocity information from the time-offlight counters, ACC, and CDC (dE/dx) [9]. The kaon identification efficiency is above 80%, while the probability of misidentifying a pion as a kaon is below 10%.
Low-momentum charged tracks that curl up in the CDC can be reconstructed multiple times by the track finder. To ensure that no track is included more than once, criteria similar to those of Ref. 10 are used. 4 In studying a mode that has a J/ψ or ψ ′ in the final state, a key strategy is to reconstruct the J/ψ only in its decays to e + e − or µ + µ − , and the ψ ′ only in its decays to e + e − , µ + µ − , or J/ψπ + π − . Although this choice abandons all but 12% of J/ψ's and 5% of ψ ′ 's, it reduces continuum backgrounds to a negligible level. The decays J/ψ → µ + µ − and ψ ′ → µ + µ − are reconstructed by combining oppositely-charged muon candidates. The invariant-mass distribution is then fitted, modeling the J/ψ and ψ ′ as double Gaussians. Muon pairs are discarded unless they have an invariant mass within ±3σ of the fitted J/ψ or ψ ′ mean, where σ is the width of the narrower Gaussian. Similarly, J/ψ → e + e − and ψ ′ → e + e − decays are reconstructed by combining oppositelycharged electron candidates. To account for energy losses due to final-state radiation or bremsstrahlung in the detector, any photons detected within 50 mrad of the initial direction of an electron candidate are also included in the e + e − invariant-mass calculation. Electron pairs are discarded unless they have an invariant mass within the range extending from −4σ to +3σ of the fitted J/ψ or ψ ′ mean. This mass window is asymmetric about the mean so as to include the radiative tails of the J/ψ and 3 Lists of four-vectors for a given decay chain are generated using EvtGen [7]. The detector response is then simulated using GEANT [8], combining randomly-triggered data with the simulated events. 4 See Ref. 11 for a detailed description of the event selection. ψ ′ , which are not completely recovered by the photon addition.
Lepton track pairs that survive the mass requirements are fitted to a common vertex, which is constrained within errors to the measured interaction point. This vertex is then fixed, and another fit is performed, constraining the dilepton invariant mass to the nominal J/ψ or ψ ′ mass. Since the observed widths of the J/ψ and ψ ′ are dominated by measurement error, this procedure improves the mass resolution of the B candidate.
To reconstruct B-meson candidates, each J/ψ or ψ ′ candidate is combined with a kaon candidate and two oppositely-charged pion candidates. Kaon and pion candidates are charged tracks that satisfy identification criteria for kaons and pions, respectively, and have an impact parameter with respect to the fitted dilepton vertex of |dr| < 0.4 cm and |dz| < 1.5 cm. 5 Any pion candidate that is identified as the product of a K 0 S → π + π − decay is discarded. 6

A. B-Meson reconstruction
Two kinematic variables can be used to identify B mesons. First, the reconstructed mass of a true B meson is likely to fall near the nominal B mass. Second, as B mesons are produced in the reaction the energy of each B in the Υ(4S) frame is half the total energy of the electron and positron beams in this frame. Since beam-energy drifts can cause the mass of the Υ(4S) to vary, it is customary to recast these kinematic variables in forms that are readily corrected for drifts in the beam energy-namely, the energy difference ∆E and beam-constrained mass M bc : 5 The impact-parameter requirement is not applied to the pions in ψ ′ → J/ψπ + π − . 6 To reconstruct K 0 S → π + π − , oppositely-charged pion candidates are combined, and the K 0 S selection criteria of Ref. 13 are applied. Both pions are vetoed if their combined invariant mass lies between 0.482 GeV/c 2 and 0.510 GeV/c 2 , which corresponds roughly to a region extending from −4σ to +3σ around the nominal K 0 S mass.
Here, P * (B) is the momentum of the B candidate in the Υ(4S) frame, while E * beam is half the energy of the Υ(4S) and is measured independently. For a correctlyreconstructed B meson, ∆E peaks at zero, and M bc peaks at the nominal B mass.
In the case of a multiparticle final state such as J/ψK + π + π − or ψ ′ K + π + π − , multiple B candidates can pose a challenge. If a correctly-reconstructed B candidate includes a low-momentum pion, then an additional B candidate can be formed by replacing that pion with a low-momentum pion from the other B. As the exchange does not significantly affect the energy or momentum of the B candidate, both candidates can satisfy ∆E and M bc criteria. Multiple candidates can spoil branching-fraction measurements and distort observed mass spectra. To ensure that no B decay is counted more than once, a best-candidate selection is performed. First, B candidates are required to have |∆E| < 0.2 GeV and M bc > 5.27 GeV/c 2 . This leaves 25% of B + → J/ψK + π + π − events and 34% of B + → ψ ′ K + π + π − events with multiple candidates; these events have a mean multiplicity of 2.4 and 2.7, respectively. If a given event has multiple B candidates with the same final state, the charged tracks that make up each B candidate are fitted to a common vertex. The candidate whose vertex fit has the smallest χ 2 is selected. According to MC studies, this procedure identifies the correct B candidate in approximately 55% of cases where there are multiple candidates.

B. Signal and sideband regions
Data distributions of ∆E for B + → J/ψK + π + π − and B + → ψ ′ K + π + π − are shown in Fig. 1. The signal is modeled as a double Gaussian with a single mean, fixing the width and relative height of the wider Gaussian to the results of a MC fit. The background is modeled as a first-order polynomial. Based on these fits, the signal region is defined as where µ ∆E is the mean of the signal peak, and σ ∆E is the width of the narrower Gaussian. The sideband region, which is used to estimate the background under the signal, is defined as 7 The small contribution from B + → X(3872)K + is not vetoed.
The sideband normalization factor is given by where p bkg is the polynomial representing the background. The fraction of signal-region events that are background is estimated as where S and B are the numbers of events in the signal and sideband regions, respectively. FIG. 1. Data ∆E distributions for B + → J/ψK + π + π − (top) and B + → ψ ′ K + π + π − (bottom). The curves show the results of the fits described in the text. Dashed and dotted lines indicate the signal and sideband regions, respectively.

IV. COORDINATE TRANSFORMATIONS
The data in the sideband region are used to model the background in the signal region. Figure 2, which shows the distribution of M (Kππ) for B + → J/ψK + π + π − events in the signal and sideband regions, reveals a problem: signal and sideband data have different end points in M (Kππ).
Plotting ∆E versus M (Kππ) reveals the cause of the discrepancy. As Fig. 3 demonstrates, the kinematically allowed range of M (Kππ) depends on ∆E. While the minimum value of M (Kππ) is M (Kππ) min = M K +2M π , the maximum value, which is attained when both the K + π + π − system and the J/ψ are at rest in the Bcandidate's rest frame, varies with ∆E as M (Kππ) max = ∆E + M B − M J/ψ . Here, M K , M π , M B and M J/ψ stand for the nominal masses of the subscripted particles. FIG. 2. Invariant mass of the K + π + π − system in B + → J/ψK + π + π − data. Open and filled histograms show events in the signal and normalized sideband regions, respectively. Transforming M (Kππ) as follows removes its dependence on ∆E: Here, Figure 4    Transformed invariant mass of the K + π + π − system in B + → J/ψK + π + π − data. Open and filled histograms show events in the signal and normalized sideband regions, respectively.
An important feature of the transformation is that it does not change M (Kππ) at ∆E = 0. Thus, although sideband and signal regions are both transformed, the change is minimal in the signal region.
Just as the range of M (Kππ) depends on ∆E, the ranges of M (Kπ) and M (ππ) also depend on ∆E. To correct for this dependence, transformations similar to Eq. 10 are applied. The variable M (Kπ) is transformed as and the variable M (ππ) is transformed as Here, M (Kπ) min = M K + M π , and M (ππ) min = 2M π .  The lines are defined as in Fig. 3. The concentration of events near 0.9 GeV/c 2 represents random particle combinations containing real K * (892)s.
As Fig. 7 illustrates, transforming the M (Kπ) coordinate distorts the shapes of the K * (892) and D 0 backgrounds. This distortion must be taken into account in parametrizing the background for the three-dimensional fits of Sec. VI (i.e., in Eqs. 24 and 25). Modeling the distortion is straightforward. First, the peak is described as a Breit-Wigner or Gaussian in the untransformed coordinate, M (Kπ). Using Eq. 11, M (Kπ) is then written as a function of M ′ (Kπ) and ∆E. The expression is numerically integrated over the relevant range of ∆E to obtain the shape of the peak as a function of M ′ (Kπ). Figure 8 demonstrates the transformation of the K * (892) background shape.
The data also contain K 0 S and ρ 0 backgrounds, albeit less prominently. The distortion of these peaks by the M (ππ) transformation is modeled by describing the K 0   the help of Eq. 12, and integrating this over the appropriate region of ∆E.
As the main source of background in this analysis is misreconstructed B-meson decays, the transformations were checked by analyzing a generic-MC simulation of Υ(4S) decays to B + B − and B 0B0 , with all known decay modes included. The M (Kππ), M (Kπ), and M (ππ) distributions were found to display the same ∆E-dependence in MC simulation as in data. Excluding signal events from the MC sample, the distributions of background events in the signal and sideband regions were compared with and without the transformations. The transformed sidebands were found to reproduce the shape of the background in the signal region more accurately than the untransformed sidebands, especially at high M (Kππ), M (Kπ), and M (ππ). 9 Although some discrepancy was observed between the background in the signal and sideband regions near the K 0 S and ρ masses, this is mostly independent of the transformation and is taken into consideration when calculating systematic errors.
The transformations of Eqs. 10-12 were also applied to B + → ψ ′ K + π + π − , with M J/ψ replaced with M ψ ′ . The results of the checks were the same.

V. TOTAL BRANCHING FRACTIONS
Branching fractions for B-meson decays to J/ψK + π + π − and ψ ′ K + π + π − final states are measured using a background-subtraction technique. 10 For each final state, data events in the signal and sideband regions are distributed into cubic bins in M 2 (Kππ), M 2 (Kπ), and M 2 (ππ). The number of signal events observed in each bin is calculated as where S i and B i are the numbers of signal-region and sideband-region events, respectively, that fall into the ith bin, and f B is the sideband normalization factor given by Eq. 8. The fraction of charged B mesons that decay to the final state in question can be expressed as where ε i is the signal efficiency in bin i, and N B is the total number of charged B mesons in the data sample. Assuming equal rates for Υ(4S) → B + B − and Υ(4S) → B 0B0 , N B is equal to the number of B pairs produced, which is measured independently.
To determine the signal efficiency, we generate 10.7 × 10 6 nonresonant B ± decays to each of the two final states of interest. We then reconstruct these signal-MC events, applying the same event-selection requirements as with data. We bin the generated events according to the generated values of M 2 (Kππ), M 2 (Kπ), and M 2 (ππ), and the reconstructed events according to the reconstructed values of M 2 (Kππ), M 2 (Kπ), and M 2 (ππ). The efficiency in each bin is the ratio of reconstructed to generated events in that bin. Figure 9 shows the dependence of the efficiency on the three variables. Figure 10 shows the corresponding data distributions. The overall efficiency is (19.85 ± 0.01)% 11 for B + → J/ψK + π + π − and (6.58 ± 0.01)% for B + → ψ ′ K + π + π − . The number of efficiency-corrected signal events observed is (4.14 ± 0.06)×10 4 for B + → J/ψK + π + π − and (1.12±0.05)×10 4 for B + → ψ ′ K + π + π − . This method of measuring branching fractions automatically corrects for efficiency variations over the phase space. It also makes no assumptions as to the shape of the signal in ∆E.

A. Systematic errors
The systematic error in the branching fractions is estimated by adding in quadrature various contributions, which are assumed to be uncorrelated. Where possible, a correction is applied.
Since we use MC simulation to determine the signal efficiency in Eq. 14, any discrepancy in signal-reconstruction efficiency between data and simulation will result in a systematic error. Based on studies of the trackreconstruction efficiency, we include a systematic error of 1.0% for each lepton track, 1.4% for each pion track, and 1.2% for each kaon track, adding linearly. Based on studies of the lepton-identification efficiency, which show that the simulation underestimates the lepton-identification efficiency, we apply a correction factor of 0.984±0.019 for each electron track and 0.962±0.031 for each muon track. Based on studies of the kaon identification efficiency, we include a systematic error of 1% for B + → J/ψK + π + π − and 2% for B + → ψ ′ K + π + π − .
The bin size of 0.15 GeV 2 /c 4 is chosen based on the dependence of the efficiency-corrected signal yield on the bin size. The error associated with this choice is taken to be the rms of the signal yield in the region between 0.1 GeV 2 /c 4 and 0.2 GeV 2 /c 4 .
In the ∆E distributions for B + → J/ψK + π + π − and B + → ψ ′ K + π + π − nonresonant MC simulation, shown in Fig. 11, a small polynomial background can be seen under the peak. Since all the events in the MC sample include a signal decay, this "background" is made up of misreconstructed signal events. Although these events are included as signal in the efficiency calculation, they are removed by the background-subtraction procedure. The fraction of signal that is subtracted in this way is found to be (3.75 ± 0.91)% for B + → J/ψK + π + π − , and (5.7 ± 1.2)% for B + → ψ ′ K + π + π − . The observed  branching fractions are corrected for this effect, and the associated uncertainty is included as a systematic error.
To determine the sideband normalization factor f B in Eq. 13, the data ∆E distribution is fitted as described in Sec. III B. In this fit, the background under the signal is parametrized as a first-order polynomial. To estimate the error introduced by this assumption, a second fit is performed, parametrizing the background as a secondorder polynomial. The fractional change in the signal yield is taken as a systematic error.
Since the signal and sideband regions are defined based on the results of fitting the data ∆E distribution, µ ∆E and σ ∆E in Eqs. 6 and 7 are varied within the fit errors.
In the MC sample used for determining the efficiency, J/ψ's from the signal B are forced to decay to e + e − or µ + µ − , and ψ ′ 's from the signal B are forced to decay to e + e − , µ + µ − , or J/ψπ + π − . Thus, to obtain branching fractions for B + → J/ψK + π + π − and B + → ψ ′ K + π + π − , the branching fractions measured using Eq. 14 are divided by previously-measured values [14] of these J/ψ and ψ ′ decay rates. The uncertainties of these previous measurements are included as a systematic error.
Finally, the error in N B is 1.3%. Table I lists the components of the systematic errors.

B. Results
The measured branching fractions are As a cross-check, we also measure a branching fraction for B + → ψ ′ K + , using a similar method but reversing the ψ ′ → J/ψπ + π − veto in the reconstruction of B + → J/ψK + π + π − . This branching fraction is which is consistent with the previously-measured value of (6.48 ± 0.35) × 10 −4 [14]. Our B + → J/ψK + π + π − branching-fraction measurement represents a significant improvement over previous measurements [14]. It is consistent with Ref. 15 but inconsistent with Ref. 16 at the 3.4-σ level. Our B + → ψ ′ K + π + π − branching-fraction measurement is also a significant improvement over the previous measurement [17].

A. Fitting technique
Signal-region data are fitted by maximizing 12 the loglikelihood function, which is given by where the sum is over the events in the signal region, x i is the vector of coordinates for a given event (i.e., , a is the vector of parameters with respect to which ℓ is maximized, and p is the probability-density function (PDF) that is used to model the observed distribution. The distribution of events in the signal region is modeled as where p B and p S describe the observed shapes of the background and signal, respectively. The constants n B and n S are the background and signal fractions in the signal region; the former is given by Eq. 9, and the latter is 1 − n B . 13 The observed signal distribution p S is expressed as 12 Standalone MINUIT [18] is used for all maximizations in this section. 13 The background fraction n B is corrected for the oversubtraction effect described in Sec. V A.  where ε is the detector efficiency, φ is the phase-space density, and s is the raw signal function.
Using nonresonant MC simulation, we have measured the detector resolution to be approximately 3-4 MeV/c 2 in each of the three coordinates M (Kππ), M (Kπ), and M (ππ). Since this is smaller than the width of any resonance included in the fits, we neglect the effect of detector resolution on line shapes.
The following five sections describe the methods followed in performing the integrals of Eq. 16 and in obtaining the functions p B ( x), ε( x), φ( x), and s( x; a) in Eqs. 16 and 17.

B. Normalization procedure
The integrations of Eq. 16 are performed numerically, using Simpson's rule. A step size of 0.010 GeV 2 /c 4 for B + → J/ψK + π + π − and 0.005 GeV 2 /c 4 for B + → ψ ′ K + π + π − is used in each dimension. 14 The threedimensional region of integration can be determined by noting that the minimum and maximum values of M 2 (Kππ) are given by where M B , M K , M π , and M ψ are the nominal values of the subscripted particles. For a given value of M 2 (Kππ), the minimum and maximum values of M 2 (ππ) are 14 The larger step size is necessary for B + → J/ψK + π + π − because of the larger phase space, which significantly increases the CPU time required for the integration.  Figure 12 shows the calculated kinematic boundaries for B + → J/ψK + π + π − , along with the observed distributions of sideband data, for a slice in M 2 (Kππ). Events that do not fall within the calculated boundaries are excluded from the fits. Because of the coordinate transformations of Sec. IV, such events are rare: 36 of the 12913 sideband events and 3 of the 10594 signalregion events for B + → J/ψK + π + π − , and 3 of the 2230 sideband events and none of the 1176 signal-region events for B + → ψ ′ K + π + π − fall outside the boundaries.

C. Background functions
To determine the three-dimensional shape of the background in the signal region (i.e., p B ( x) in Eq. 16), an unbinned maximum-likelihood fit is performed on the sideband-region data. The log-likelihood function to maximize is given in this case by where the sum is over the events in the sideband region. The maximization is performed by varying the parameters a B , which are then fixed at their optimal values in fitting the signal region. The background is modeled as the sum of a combinatorial term and a set of noninterfering resonances. For and for B + → ψ ′ K + π + π − , In Eqs. 24 and 25, T n represents an nth-order Chebyshev polynomial. The variables x, y, and z stand for M 2 (Kππ), M 2 (Kπ), and M 2 (ππ), respectively, and are defined over the intervals The peak functions BW K * (892) , G D , G KS , and BW ρ are obtained as described in Sec. IV. Each peak function P ( x) is normalized over the kinematically-allowed phase space to satisfy The factor of e −2x that modulates the peak functions was found empirically to produce a good fit to the sideband data. 15 Table II lists the fitted parameters of the background functions. The statistical error in each parameter is defined as the change in that parameter required to reduce the log likelihood by 1/2. The fitted functions, normalized to the total number of events in the fit, are shown projected onto the three axes along with the sideband data in Fig. 13. Figures 14 and 15 show M 2 (Kπ) and M 2 (ππ) projections for slices in M 2 (Kππ).
As a measure of goodness of fit, a χ 2 variable is calculated by distributing the data into cubic bins that are 0.1 GeV/c 2 wide on each side. The normalized PDF, with the parameters set to their best-fit values, is integrated over each bin and multiplied by the total number of events in the fit to determine the number of events expected in the bin. Adjacent bins are combined until each bin has at least 6 data events. A χ 2 variable for the multinomial distribution is then calculated as [19] where N bins is the total number of bins used, n i is the number of observed events in a given bin, and p i is the number expected in that bin based on the PDF. If the expected distribution p i were obtained by a binned maximum-likelihood fit of the data distribution n i , the number of degrees of freedom associated with this χ 2 would be reduced by the number of fit parameters Results of sideband fits for B + → J/ψK + π + π − (left) and B + → ψ ′ K + π + π − (right). Data (points) and fits (histograms) are shown projected onto the three axes. The red histograms show the overall background functions. The combinatorial components are shown in gray, while the K * (892), ρ, KS and D backgrounds are shown in blue, magenta, cyan, and green, respectively. The K * (892) and ρ peaks are broader in B + → ψ ′ K + π + π − than in B + → J/ψK + π + π − because the distortion shown in Fig. 8 is larger in the former mode.
N par and would be given by N DOF = N bins − N par − 1.
If, on the other hand, the two distributions were not correlated by a fit, the number of degrees of freedom would be N DOF = N bins − 1. Since, in this case, the distributions are related by an unbinned maximum-likelihood fit, the true N DOF can be expected to lie between these extremes [20]. For the B + → J/ψK + π + π − sideband-data fit, ℓ B = −21484.6, while χ 2 = 1709.5 with N bins = 1707 and N par = 12. For the B + → ψ ′ K + π + π − sideband-data fit, ℓ B = 822.7, while χ 2 = 286.8 with N bins = 294 and N par = 5.

D. Efficiency functions
The dependence of the detector efficiency on the kinematic variables (i.e., ε( x) in Eq. 17) is obtained for threedimensional bins, 0.15 GeV 2 /c 4 -wide on each side, using nonresonant signal-MC simulation as described in Sec. V and illustrated in Fig. 9. The function is implemented as a lookup table: the efficiency for a given data point is the efficiency in the corresponding bin.

E. Phase-space densities
Four-body phase-space densities (i.e., φ( x) in Eq. 17) for B + → J/ψK + π + π − and B + → ψ ′ K + π + π − are obtained by using GENBOD [21] to generate final-stateparticle four-momenta that are weighted by the density of states in phase space [22]. For each decay mode, 10 8 events are generated. Event phase-space weights are distributed into cubic bins in M 2 (Kππ), M 2 (Kπ), and M 2 (ππ), with a bin width of 0.02 GeV 2 /c 4 . The phase-space density is implemented as a lookup table: the value of φ( x) for a given data point is the total phasespace weight in the corresponding bin. 16 In Fig. 16, the three-dimensional histogram of phase-space weights is projected onto the three axes, showing the distribution that signal events would have in the absence of resonant effects. Figure 16 does not indicate the functional form of φ( x), since the projection onto a single dimension effectively in-   FIG. 15. B + → ψ ′ K + π + π − sideband data (points) and fit results (histograms) for slices in M 2 (Kππ). The fit components are color coded as in Fig. 13. tegrates over the other two dimensions, and the region of integration is the complicated one described in Sec. VI B. In Fig. 17, the same projections are performed over a narrow slice in each of the other two dimensions, to illustrate the dependence of the function φ( x) on each variable.

F. Signal functions
The K + π + π − final state is modeled as a nonresonant signal plus a superposition of initial-state resonances R 1 . The latter are assumed to decay through intermediatestate resonances R 2 as R 1 → aR 2 , R 2 → bc, where a, b, and c are the final-state particles. Specifically, the function s( x; a) of Eq. 17 is expressed as Here, J 1 and J 2 stand for the spin-parity (J P ) of R 1 and R 2 , respectively. Resonances with different J 1 are added incoherently, while those with the same J 1 are added coherently. The parameters varied in the fit are the complex coefficients a nr and a J1J2 , collectively referred to as a k . While the nonresonant signal is assumed to be constant over the phase space, the resonant decay amplitudes A J1J2 are expressed as where Γ R2 (m bc ) is the mass-dependent width and F R is the Blatt-Weisskopf barrier factor The meson radial parameter R is set to 1.5 (GeV/c) −1 .
The function α J1J2 describes the spin-dependent angular distribution of the final state and is shown for various combinations of J 1 and J 2 in Table III. Resonances with spin greater than two are not included in the fitting model. In cases where there is more than one covariant spin amplitude, only the lowest spin is included. In Eqs. 30-32 and Table III, the nominal masses of the resonances R 1 and R 2 are denoted by M R1 and M R2 , and the nominal widths by Γ R1 and Γ R2 . The angle θ is between a and b in the bc rest frame and can be expressed as The variable z is given by The breakup momentum p is the momentum of a or bc in the abc rest frame: while q is the momentum of b or c in the bc rest frame: where the constant q 0 is the value of q evaluated at m bc = M R2 .
Since the components of the signal function are not individually normalized, it is not meaningful to compare  , and M 2 (ππ), for B + → J/ψK + π + π − (left) and B + → ψ ′ K + π + π − (right). In each case, the region selected is indicated above the plot. the moduli of the complex coefficients a k in Eq. 28. A decay fraction is therefore calculated for each component by integrating the component over the kinematicallyallowed region and dividing by the integral of the full signal function The integrations in Eq. 37 are performed as described in Sec. VI B. Because of interference effects, decay fractions for a given final state will not, in general, add up to unity.

G. Statistical errors
As with the sideband-region fits, the statistical uncertainties in the fit parameters (i.e., moduli and phases) are determined by the fitter: the error in a given parameter is the change in that parameter that reduces the log likelihood by 1/2. The statistical uncertainties in the decay fractions, on the other hand, are more complicated. Since a given decay fraction involves the integral of the full signal function, the error in a single decay fraction incorporates the errors in all of the parameters. To determine the statistical errors in the decay fractions, 1000 sets of correlated signal-function parameters are drawn from Gaussian distributions using the fitted parameter values and the error matrix. 17 Decay fractions are calculated for each set of generated parameters. The rms of the resulting distribution provides an estimate of the statistical error in the decay fraction.

H. Systematic errors
Several sources of systematic error are considered, as described below. They are added in quadrature to obtain the systematic errors reported in Sec. VI I.

Background parametrization
A possible source of systematic error in the fits is the fixed background fraction n B in Eq. 16. While the error in n B is small, the correction for the oversubtraction, described in Sec. V A, lowers n B by 10.8% for B + → J/ψK + π + π − and by 11.4% for B + → ψ ′ K + π + π − . The systematic error associated with this correction is estimated conservatively as the change in each parameter when the fits are performed with the uncorrected values of n B .
There may be an additional systematic error if the background in the signal region is not correctly 17 Correlated Gaussian distributions are generated using CORSET and CORGEN [24]. parametrized by the shape determined by fitting the sidebands. As noted in Sec. IV, generic-MC studies suggest that not enough of the K 0 S and ρ background peaks are removed by the sideband subtraction. To estimate this error, a fit is performed in which the coefficients of the background peaks in Eqs. 24 and 25 are doubled.

Efficiency
To estimate the error introduced by binning the efficiency information, the fits are repeated using bin sizes of 0.10 GeV 2 /c 4 and 0.20 GeV 2 /c 4 for the efficiency. The average absolute change in each parameter is the estimate of the error.
Another possible source of error is that the MC simulation may not faithfully reproduce the detector efficiency for low-momentum particles. To test for such an effect, two additional fits are performed. In the first fit, only charged particles with a momentum greater than 200 MeV/c are included. In the second fit, the |dr| and |dz| requirements described in Sec. III are loosened from 0.4 cm to 0.8 cm, and from 1.5 cm to 3.0 cm, respectively. The changes in each parameter observed in these two fits are added in quadrature to obtain an estimate of the error due to inaccuracies in the efficiency estimation.
Using only the three variables M 2 (Kππ), M 2 (Kπ), and M 2 (ππ) in this analysis is equivalent to integrating over variables that describe the relative momentum of the J/ψ or ψ ′ with respect to the K + π + π − system. In this integration, the terms corresponding to K + π + π − states with different initial-state spin-parity cancel out, producing Eq. 28. This cancellation, however, is exact only if the detector efficiency is flat over the extra variables. To determine the effect of neglecting these variables, an additional set of fits is performed, in which the efficiency in Eq. 17 is calculated as a function of the two angles between the J/ψ or ψ ′ and the K + π + π − system, rather than M 2 (Kππ), M 2 (Kπ), and M 2 (ππ). The resulting fitted parameters are compared to those obtained by a fit in which the efficiency is held constant. 18 The absolute change in each parameter is found to be small (less than 15% of the statistical error) and is included in the systematic error.

Integration step size
To estimate the error introduced by the finite step size used in the numerical integrals of Secs. IV and VI B, the fits are repeated, using a step size of 0.005 GeV 2 /c 4 for B + → J/ψK + π + π − and 0.010 GeV 2 /c 4 for B + → ψ ′ K + π + π − . The change in each parameter is an estimate of the uncertainty associated with the numerical integration.

Modeling of the signal
The masses and widths of the resonances included in the fits are listed in Table IV. To estimate the systematic error associated with the uncertainties in these quantities, the fits are repeated, varying each fixed quantity within its errors. For each mass or width, the average absolute change in each parameter is recorded. These average changes are then added in quadrature.
In fitting the B + → J/ψK + π + π − data, the modulus for K * 2 (1430) → K * (892)π is allowed to float. Relative to this modulus, the moduli 19 for K * 2 (1430) → Kρ and K * 2 (1430) → Kω are fixed based on previously-measured relative branching fractions [14]. To estimate the associated systematic error, additional fits are performed, varying these branching fractions within their uncertainties. Table V lists the values of the moduli and phases of the complex coefficients a k of Eq. 28 obtained by fitting signal-region data for B + → J/ψK + π + π − , as well as the corresponding values of the decay fractions, given 19 The phases of the three submodes are allowed to float. by Eq. 37. The fitted PDF is shown projected onto the three axes, along with the data, in Fig. 18. Figure 22 shows M 2 (Kπ) and M 2 (ππ) projections for slices in M 2 (Kππ). The legend is presented in Fig. 19. In this fit, ℓ = −10575.4, while χ 2 = 1475.9 with N bins = 1202 and N par = 28. Similarly, Table VI shows the fitted parameters for B + → ψ ′ K + π + π − signal-region data, as well as the corresponding decay fractions. Figure 20 shows the fitted PDF and data projected onto the three axes, while Fig. 23 shows M 2 (Kπ) and M 2 (ππ) projections for slices in M 2 (Kππ). In this fit, ℓ = 638.3, while χ 2 = 180.1 with N bins = 168 and N par = 10.

I. Results
Finally, the B + → J/ψK + π + π − signal-region data are fitted again, this time floating the mass and width of the (39) Table VII shows the fitted parameters, along with the corresponding decay fractions. Figure 21 shows the fitted PDF and data projected onto the three axes, while    tematic error, which is somewhat offset by an increase in the corresponding statistical error. In particular, the K 1 (1270) → K * (892)π decay fraction is especially sensitive to the K 1 (1270) mass and width. In any fit involving many floating parameters, local likelihood maxima can present a problem. To ensure that the fit results are global maxima, 100 additional fits were performed for each of the three cases, selecting random starting values for the parameters. None of these fits yielded better likelihoods than those presented above. The local maxima encountered in the course of this test are discussed in the Appendix.
J. Discussion

Signal components
In choosing the signal components to be included in the fits, the data were used as a guide. As the K 1 (1270) signal is prominent in both B + → J/ψK + π + π − and TABLE V. Fitted parameters of the signal function for B + → J/ψK + π + π − , along with the corresponding decay fractions.
As a further guide, the decays B 0 → J/ψK + π − and B 0 → ψ ′ K + π − were reconstructed. The observed Kπ mass spectra are shown in Fig. 25. Consistent with the 1 + spin-parity assignment of the K 1 (1270), no K 1 (1270) → Kπ signal appears in these spectra. In both modes, a small peak can be seen near 1.4 GeV/c 2 in M (Kπ); this may have contributions from K * (1410) or K * 2 (1430), as well as K * 0 (1430). The absence of a K * (1680) peak in B 0 → J/ψK + π − is noteworthy, although a precise statement would require an analysis of the efficiency and phase space for these modes. 20 In B 0 → ψ ′ K + π − , the kinematically-allowed M (Kπ) region does not allow any conclusions to be drawn about the presence or absence of a low K * (1680) tail. Observed Kπ mass spectra for B 0 → J/ψK + π − (top) and B 0 → ψ ′ K + π − (bottom) data.

Interference effects
The inclusion of interference among submodes sharing the same initial-state spin-parity is essential to obtaining good fits to the data. In particular, dramatic interference effects are observed between K 1 (1270) → K * (892)π and K 1 (1270) → Kρ, as well as between K 1 (1270) → Kρ and K 1 (1270) → Kω. 20 A detailed analysis of B → J/ψKπ and B → ψ ′ Kπ is beyond the scope of this work. A Dalitz analysis of the latter mode was presented in Ref. 25. Figure 27 shows scatterplots of signal-region B + → J/ψK + π + π − data over the three coordinates. Interference between K 1 (1270) → K * (892)π and K 1 (1270) → Kρ is responsible for the weakening of the latter signal at M (Kπ) > M K * (892) . Although the four-body phase space decreases with increasing M (Kπ), this is not sufficient to account for the abrupt falloff. To describe the data in this region, the two modes must be added coherently.
Since the previously-measured [14] branching fraction for K 1 (1270) → Kω is small compared to that for K 1 (1270) → Kρ, and since only 1.5% of ω's decay to π + π − , 21 one might expect K 1 (1270) → Kω to play a negligible role in this analysis. Nonetheless, since the ω is much narrower than the ρ, it significantly distorts the observed ρ line shape through interference [27]. In Fig. 26, the M 2 (ππ) projections of Figs. 18, 20 and 21 are finely binned to demonstrate this interference pattern, which is accurately modeled by the PDFs. Although ω decays dominantly to to π + π − π 0 , it can also decay to π + π − through G-parity violation [26], which causes mixing between ρ and ω. An ω component is therefore present whenever a particle decays to ρ.
The peculiar shape of the observed ρ-ω interference pattern is caused by kinematic effects. The largest contribution to the ρ signal comes from K 1 (1270) → Kρ, which straddles the edge of phase space, as can be seen in the middle panel of Fig. 27. The distortion that is caused by this kinematic cutoff is taken into account automatically by integrating the signal function only over the kinematically-allowed region, as described in Sec. VI A. Modeling the data accurately requires including ρ-ω interference, incorporating the four-body phase space factor into the signal function, and integrating the signal function over only the kinematically-allowed phase space.

The L region
The M (Kππ) region between 1.5 and 2.0 GeV/c 2 , historically referred to as the L region, comprises several wide, overlapping resonances [28][29][30]. The large uncertainties in the masses and widths of the known states in this region make it difficult to characterize this region in this analysis. The model presented here is not necessarily the only one supported by the data.
To describe the structure observed at 2.6 GeV 2 /c 4 in the M 2 (Kππ) distribution of B + → J/ψK + π + π − , a peak with a mass of 1.605 GeV/c 2 and a width of 115 MeV/c 2 is included in the fit, decaying to K * (892)π and Kρ. This peak, which is referred to as K(1600) in this paper, may be the K 2 (1580), an as-yet unconfirmed J P = 2 − state that has previously been observed decaying to K * (892)π [29].
As can be seen in Fig. 22, the high end of the M 2 (Kππ) spectrum of B + → J/ψK + π + π − exhibits K * and ρ signals. To fit the data in this region, we include a K * 2 (1980) resonance, which is another state that currently requires confirmation.
Fitting the B + → ψ ′ K + π + π − data is more difficult still, as there are fewer events to analyze, and only a small portion of the L-region is within the kinematic limits of the decay. In addition to the K 1 (1270) signal, the M 2 (Kππ) spectrum contains what appears to be the low-mass tail of at least one high-mass resonance. As Fig. 23 shows, there are clear K * (892) and ρ peaks at high M 2 (Kππ); these are not reproduced by the PDF if no high-mass resonance is included in the model. If the enhancement is modeled as a single resonance, the data favor a mass of roughly 1.7 GeV/c 2 and a width of 400-500 MeV/c 2 . In this analysis, the enhancement is modeled as the K * (1680). The data do not preclude other possibilities, such as the K 2 (1770). Indeed, the hint of f 0 (980) in the last M 2 (Kππ) slice in Fig. 23 cannot come from a 1 − state such as the K * (1680), or from a 2 + state such as the K * 2 (1430).

Comparison with previous measurements
It is interesting to compare the relative decay fractions for K 1 (1270) submodes in the B + → J/ψK + π + π − fits to previous measurements of K 1 (1270) branching fractions. For this purpose, we use the decay fractions with phase space shown in Tables V and VII, include isospin factors, and assume branching fractions of (1.53 +0. 11 −0.13 )% for ω → π + π − , and (93±10)% for K * 0 (1430) → Kπ [14]. The calculation neglects the systematic errors in the decay fractions and assumes that the statistical errors among the decay fractions are uncorrelated. Moreover, it assumes that the K 1 (1270) decays only to K * (892)π, Kρ, Kω, and K * 0 (1430)π, and neglects interference among these decay channels. The comparison is shown in Table VIII. While the ratios of the K 1 (1270) branching fractions to K * (892)π, Kρ, and Kω are consistent with the previously-measured values, the branching fraction to K * 0 (1430)π is significantly smaller. As shown in Sec. VI I, the data favor a smaller mass and a larger width for the K 1 (1270) than the Particle Data Group (PDG) values. This is mainly due to the excess of K * (892) and ρ at low M 2 (Kππ), as can be ascertained by comparing the first row of plots in Figs  . The effect is most apparent in the left and right plots. Scatterplots for B + → ψ ′ K + π + π − data are similar but limited by statistics.

Limitations of the method
There are large uncertainties in the masses and widths of many of the states included in the fits, as can be seen in Table IV. Although this is taken into account in calculating the systematic error, it nonetheless limits the accuracy of the model.
In B + → ψ ′ K + π + π − , the small sample size and the kinematic cutoff limit the conclusions that can be drawn about the signal components. In B + → J/ψK + π + π − , the sample size is larger, but a further limitation is imposed by the increase in computation time as more parameters are added to the fit. Each additional decay channel that is included in the signal function contributes a modulus and possibly a phase to be varied in the fit. Since the normalization integral of the signal function in Eq. 16 depends on the values of the parameters a, the integration must be performed for each set of parameters attempted by the fitter. While the step size used in the numerical integration can be increased to speed up the process, it must be small enough to allow the PDF to resolve the structures in the data. In particular, the ρ-ω interference pattern can be fitted with a step size of 0.01 GeV 2 /c 4 , but not with a step size of 0.02 GeV 2 /c 4 . As a consequence of the finite processor speed, not every possible decay channel can be included in the fit. The model is necessarily incomplete.
The large nonresonant component seen in both B + → J/ψK + π + π − and B + → ψ ′ K + π + π − may be an indication of contributions from additional wide kaon excitations. It may also incorporate some misreconstructed resonant signal. While the nonresonant component is assumed in this analysis to be distributed according to phase space, this assumption may be inaccurate. There are currently no accepted models of nonresonant Bmeson decays.
It is difficult, in an analysis like the one presented here, to determine the significance of a given component of the signal. An improvement in the likelihood upon the addition of a new resonance to the signal function indicates only that the model is incomplete, not necessarily that the data contain the particular resonance. Furthermore, unless the model is accurate in every other way, floating the mass and width of a particle in the fit may not yield a reliable result, as the fitter may set these parameters to compensate for the model's deficiencies. This is especially important in the high-M 2 (Kππ) region, where the statistics are limited and there are large uncertainties in the masses and widths of the resonances included in the signal function. Thus, although the K(1600) component of the signal function for B + → J/ψK + π + π − greatly improves the quality of the fit, it is difficult to claim that it is a single particle, let alone measure its mass and width.

VII. CONCLUSIONS
Using data recorded by the Belle detector, we have measured branching fractions for the decays B + → J/ψK + π + π − and B + → ψ ′ K + π + π − with improved precision (see Sec. V B). We have also performed amplitude analyses in three dimensions-M 2 (Kππ), M 2 (Kπ), and M 2 (ππ)-to determine the resonant structure of the K + π + π − final state in these decays (see Sec. VI I).
We have shown that the K 1 (1270), which is the dominant component of the K + π + π − final state in B + → J/ψK + π + π − , is also prominent in B + → ψ ′ K + π + π − . The large sample available for the former decay reveals a small peak at M (Kππ) ≈ 1.4 GeV/c 2 . Our threedimensional fits represent a first attempt to determine the components of this peak.
Performing an unbinned fit in three dimensions exploits practically all of the information available in the data. While it is relatively easy to obtain a good fit in one dimension, requiring a fit that succeeds in three dimensions greatly restricts the class of successful models.
With high statistics, it is possible to use interference effects and the spin-dependent angular distribution of the final state to distinguish overlapping resonances. In particular, we have shown that ρ-ω interference cannot be neglected in studying K 1 (1270) decays to Kπ + π − final states.
The large size of the B + → J/ψK + π + π − data sample allows us to measure the mass and width of the K 1 (1270) with improved precision (see Eqs. 38 and 39). These values differ considerably from previously-published values [14]. The analysis of these data also provides information on the relative strengths of K 1 (1270) decays to Kρ, Kω, K * (892)π, and K * 0 (1430)π final states (see Table VIII). While the results are consistent with previous measurements for the first three modes, they indicate a much smaller rate of decay to K * 0 (1430)π than previously accepted.
Although more data are required to clarify the structure of the high M 2 (K + π + π − ) region in both B + → J/ψK + π + π − and B + → ψ ′ K + π + π − , we have shown that this region contains broad resonances that decay to K * (892)π and Kρ final states.
The analysis presented in this paper demonstrates that the decay modes B + → J/ψK + π + π − and B + → ψ ′ K + π + π − can provide clean laboratories for the spectroscopy of excited kaon states. Many of these states still require confirmation or more precise mass and width measurements. As more data become available at future super-B factories, analyses similar to the one presented here can further elucidate the higher regions of the kaon spectrum.

ACKNOWLEDGMENTS
We thank the KEKB group for the excellent operation of the accelerator, the KEK cryogenics group for the efficient operation of the solenoid, and the KEK computer group and the National Institute of Informatics for valuable computing and SINET3 network support. We acknowledge support from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, the Japan Society for the Promotion of Science (JSPS), and the Tau Research in a Priority Area ("New Development of Flavor Physics"), and from JSPS for Creative Scientific Research ("Evolution of Tau-lepton Physics").

APPENDIX: LOCAL MAXIMA
This appendix summarizes the results of the localmaximum test, in which each of the three signal-region fits was repeated 100 times with randomly selected starting values for the parameters. In each case, the best likelihood obtained coincided with the solution presented in Tables V-VII. In the following, these solutions are referred to as the "global maxima." For the B + → J/ψK + π + π − fit with the K 1 (1270) mass and width fixed to their values in Table IV, the local maximum closest to the global maximum presented in Table V had a likelihood of −10608.6, which is 8.2σ away from the global maximum.
For the B + → ψ ′ K + π + π − fit, two local maxima were found: one with a likelihood of 635.4 and the other with a likelihood of 635.9; these are 2.4σ and 2.2σ away from the global maximum, respectively. The former had an unphysically large decay fraction for K 1 (1270) → Kω and was discarded. For the latter, all the parameters were within statistical error of the values presented in Table VI, with the exception of the K 1 (1270) → K * (892)π amplitude and decay fraction, which were higher by 1.4 times the statistical error.
For the B + → J/ψK + π + π − fit with the K 1 (1270) mass and width allowed to float, two local maxima were found, both with a likelihood of −10528.8, which is 2.7σ away from the global maximum. The fitted parameters and decay fractions for these local maxima are presented in Tables IX and X