$K^{*}(892)^0$ meson production in inelastic $p+p$ interactions at 40 and 80 GeV/$c$ beam momenta measured by NA61/SHINE at the CERN SPS

Measurements of $K^{*}(892)^0$ resonance production via its $K^{+}\pi^{-}$ decay mode in inelastic $p+p$ collisions at beam momenta 40 and 80 GeV/$c$ ($\sqrt{s_{NN}}=8.8$ and 12.3 GeV) are presented. The data were recorded by the NA61/SHINE hadron spectrometer at the CERN Super Proton Synchrotron. The \textit{template} method was used to extract the $K^{*}(892)^0$ signal. Transverse momentum and rapidity spectra were obtained. The mean multiplicities of $K^{*}(892)^0$ mesons were found to be $(35.1 \pm 1.3 \mathrm{(stat)} \pm 3.6 \mathrm{(sys)) \cdot 10^{-3}}$ at 40 GeV/$c$ and $(58.3 \pm 1.9 \mathrm{(stat)} \pm 4.9 \mathrm{(sys)) \cdot 10^{-3}}$ at 80 GeV/$c$. The NA61/SHINE results are compared with the EPOS1.99 and Hadron Resonance Gas models as well as with world data. The transverse mass spectra of $K^{*}(892)^0$ mesons and other particles previously reported by NA61/SHINE were fitted within the Blast-Wave model. The transverse flow velocities are close to 0.1--0.2 of the speed of light and are significantly smaller than the ones determined in heavy nucleus-nucleus interactions at the same beam momenta.


Experimental setup
The NA61/SHINE experiment [1] uses a large acceptance hadron spectrometer located in the North Area of the CERN accelerator complex. The schematic layout of the NA61/SHINE detector configuration (used for p+p data taking) is shown in Fig. 1. Only the detector components, which were used in this analysis, are described below. A more detailed description of the full detector can be found in Ref. [1].  Figure 1.: The schematic layout of the NA61/SHINE spectrometer (horizontal cut, not to scale) used for p+p data taking. The beam and trigger detector configuration is shown in the inset (see Refs. [23,24] for detailed description). The chosen coordinate system is drawn on the lower left: its origin lies in the middle of the VTPC-2, on the beam axis.
A set of scintillation and Cherenkov counters (S1, S2, V0, V1 p , V1, CEDAR, THC), as well as beam position detectors (BPDs) upstream of the spectrometer, provide timing reference, identification, and position measurements of incoming beam particles. The trigger scintillation counter S4 placed 3.7 meters downstream of the target is used to select events with collisions in the target area by the absence of a charged particle hit.
Secondary beams of positively charged hadrons are produced from 400 GeV/c protons extracted from the SPS accelerator. The primary proton beam was directed to the T2 target (located 535 m before the NA61/ SHINE production target) where it interacted. Then the produced hadrons were used to form a secondary proton beam with chosen beam momentum (here 40 GeV/c and 80 GeV/c).
For 40 GeV/c p+p data taking, two Cherenkov counters, a CEDAR [25] (CEDAR-W for 40 GeV/c), and a threshold counter (THC) were used to identify particles of the secondary hadron beam. For 80 GeV/c proton beam, only the CEDAR-N counter was used. The CEDAR counter, using a coincidence of six out of the eight photo-multipliers placed radially along the Cherenkov ring, provides positive identification of protons, while the THC, operated at a pressure lower than the proton threshold, is used in anti-coincidence in the trigger logic. A selection based on signals from the Cherenkov counters allowed to identify beam protons with a purity of about 99%. A consistent value for the purity was found by bending the beam into the TPCs (Time Projection Chambers) with the full magnetic field and using identification based on its specific ionization energy loss dE /dx [26].
The main NA61/SHINE tracking devices are four large volume Time Projection Chambers located behind the target. Two of them, the vertex TPCs (VTPC-1 and VTPC-2), are located in the magnetic fields of two super-conducting dipole magnets with a combined maximum bending power of 9 Tm corresponding to about 1.5 T and 1.1 T fields in the upstream and downstream magnets, respectively. This field configuration was used for p+p data taking at 158 GeV/c [9]. In order to optimize the acceptance of the detector, the field in both magnets was adjusted proportionally to the beam momentum. The VTPCs are filled with a mixture of argon and carbon dioxide in 90/10 proportion. Each of the VTPCs provides up to 72 points on the particle trajectory. Two large main TPCs (MTPC-L and MTPC-R) are positioned downstream of the magnets symmetrically to the beam line. The MTPCs are filled with a mixture of argon and carbon dioxide in 95/5 proportion. Particle trajectories in MTPC-L or MTPC-R are determined by the use of up to 90 points. The fifth small TPC (GAP TPC) is placed between VTPC-1 and VTPC-2 directly on the beam line. It closes the gap between the beam axis and the sensitive volumes of the other TPCs. The GAP TPC is filled with a mixture of argon and carbon dioxide in 90/10 proportion, and it provides up to 7 points on the particle trajectory. Particle identification in the TPCs is based on measurements of the specific energy loss (dE /dx) in the chamber gas.
The p+p data used in this analysis were recorded with the proton beam incident on a liquid hydrogen target (LHT), a 20.29 cm long cylinder situated upstream of the entrance window of VTPC-1.

Data sets
The results on K * (892) 0 production in inelastic p+p interactions at p beam =40 GeV/c and 80 GeV/c are based on data recorded in 2009. The numbers of events selected by the interaction trigger were 4.70M and 3.87M, respectively. Table 1 presents the numbers of events recorded with the interaction trigger and the numbers of events selected for the analysis (see Sec. 3.3). The drop in the number of events after cuts is caused mainly by BPD reconstruction inefficiencies and off-target interactions accepted by the trigger. The numbers of tracks, also given in Table 1, refer to tracks registered in accepted events only. The list of track cuts is discussed in Sec. 3

Analysis method
The detailed descriptions of NA61/SHINE calibration, track and vertex reconstruction procedures, as well as simulations used to correct the reconstructed data, are discussed in Refs. [23,24,27]. Below, only the specific analysis technique developed for the measurement of the K * (892) 0 spectra in p+p interactions is described. The analysis procedure consists of the following steps: (i) selection of events and tracks (details are given in Sec. 3.3 and 3.4), (ii) selection of K + and π − candidates based on the measurement of their ionization energy loss (dE/dx) in the gas volume of the TPCs (details are given in Sec. 3.5), (iii) preparation of invariant mass distributions of K + π − pairs (details are given in Sec. 3.6), (iv) preparation of invariant mass distributions of K + π − pairs for mixed events and Monte Carlo templates (details are given in Sec. 3.6), (v) extraction of K * (892) 0 signals and obtaining the raw numbers of K * (892) 0 (details are given in Sec. 3.6 and 3.7), (vi) application of corrections (obtained from simulations) to the raw numbers of K * (892) 0 ; they include losses of inelastic p+p interactions due to the on-line and off-line event selection as well as losses of K * (892) 0 due to track and pair selection cuts and the detector geometrical acceptance (details are given in Sec. 3.8 and 3.9).

Event selection
Inelastic p+p interactions, used in this analysis, were selected by the following criteria: (i) an interaction was recognized by the trigger logic (the detailed description can be found in Refs. [23,24]), (ii) no off-time beam particle was detected within ±1 µs around the trigger (beam) particle, (iii) the trajectory of the beam particle was measured in at least one plane of BPD-1 or BPD-2 and in both planes of the BPD-3 detector, (iv) the primary interaction vertex fit converged, (v) the z position (along the beam line) of the fitted primary p+p interaction vertex was found between -590 cm and -572 cm, where the center of the LHT was at -581 cm (the range of this cut was selected to maximize the number of good events and minimize the contamination by off-target interactions), (vi) events with a single, well-measured positively charged track with absolute momentum close to the beam momentum (p > p beam − 1 GeV/c) were rejected.
The event cuts listed above select well-measured inelastic p+p interactions. The background due to elastic interactions was removed via cuts (iv) and (vi). The contribution from off-target interactions was reduced by cut (v). The losses of inelastic p+p interactions due to the event selection procedure were corrected for using simulations (see below).
The numbers of events left after the above cuts were 1.34 × 10 6 and 1.26 × 10 6 for 40 GeV/c and 80 GeV/c, respectively.

Track selection
After adopting the event selection criteria, a set of track quality cuts were applied to individual tracks. They were used to ensure high reconstruction efficiency, proper identification of tracks, and to reduce the contamination of tracks from secondary interactions, weak decays, and off-time interactions. The tracks were selected according to the following criteria: (i) the track fit including the interaction vertex converged, (ii) the total number of reconstructed points on the track was higher than 30, (iii) the sum of the number of reconstructed points in VTPC-1 and VTPC-2 was higher than 15 or the number of reconstructed points in the GAP TPC was higher than 4, (iv) the distance between the track extrapolated to the interaction plane and the interaction point (socalled impact parameter) was smaller than 4 cm in the horizontal (bending) plane and 2 cm in the vertical (drift) plane, (v) the track total momentum (in the laboratory reference system) was p lab ≤ 35 GeV/c for 40 GeV/c beam momentum and p lab ≤ 74 GeV/c for 80 GeV/c beam momentum, (vi) the track transverse momentum (p T ) was required to be smaller than 1.5 GeV/c, (vii) dE /dx track cuts were applied to select K + and π − candidates (details are given in Sec. 3.5).
The numbers of tracks left after the above cuts were 1.53 × 10 6 and 2.13 × 10 6 for 40 GeV/c and 80 GeV/c, respectively.

Selection of kaon and pion candidates
In this analysis, charged particle identification was based on the measurement of ionization energy loss (dE / dx) in the gas volume of the TPCs. In Fig. 2 the example (for 40 GeV/c data) dE / dx values as a function of total momentum (p lab ), measured in the laboratory reference system, are shown for positively and negatively charged particles, separately. For both beam momenta (40 GeV/c and 80 GeV/c) the K + and π − candidates were selected by requiring their dE /dx values to be within (−1.2σ; +1.8σ) (for kaons) and (−2.7σ; +3.3σ) (for pions) around their empirical parametrizations of Bethe-Bloch curves (lines in Fig. 2). The quantity σ represents a typical standard deviation of a Gaussian function fitted to the dE /dx distribution of charged kaons and pions. Since only small variations of σ were observed for different total momentum and transverse momentum bins, fixed values σ = 0.044 for K + and σ = 0.052 for π − were used. The asymmetric cuts were applied to reduce the number of protons within kaon candidates and the number of kaons within pion candidates. Moreover, the upper limits for p lab were introduced (p beam − 5 GeV/c for 40 GeV/c, p beam − 6 GeV/c for 80 GeV/c) in order to eliminate the region where dE /dx calibration is less reliable (due to low statistics). 3.6. K * (892) 0 signal extraction The K * (892) 0 lifetime is about 4 fm/c [28], so this meson resonance decays essentially at the primary interaction vertex. The raw numbers of K * (892) 0 mesons are obtained by performing fits to backgroundsubtracted invariant mass spectra of K * (892) 0 decay products. The invariant mass is defined as m K + π − = 2 , where E represents the total energy and p the momentum vector of daughter particles from K * (892) 0 decay.
In this analysis, the template method (see below) was applied to extract the raw numbers of K * (892) 0 particles. Its advantages over the standard method (based on mixed events only) were described in Ref. [9]. The template method was already successfully used by NA61/SHINE in the analysis of K * (892) 0 production in p+p interactions at 158 GeV/c.
In the template method the invariant mass spectra of the data (blue data points in Figs. 3 and 4 (left)) were fitted with a function given by Eq. (1): The background is described as a sum of two contributions: T MC res and T DAT A mix . The T DAT A mix component is the combinatorial background estimated based on the mixing method (invariant mass spectra calculated for K + π − pairs originating from different events). The T MC res template (MC abbreviation stands for Monte Carlo) is the shape of the simulated background, which describes the contribution of K + π − pairs originating from (i) combination of tracks that come from decays of resonances different than K * (892) 0 , for example, one track from a ρ 0 meson and one from a K * + meson, (ii) combination of tracks where one comes from the decay of a resonance and one comes from direct production in the primary interaction.
The MC samples used to prepare the T MC res templates were generated by the Epos1.99 [29] hadronic interaction model using the CRMC 1.4 package [30]. Generated p+p events were processed through the NA61/SHINE detector simulation chain and then through the same reconstruction routines as the data. The MC simulation maintains the history of particle production, thus allowing to check their identity and origin, enabling the construction of the proper templates. For the reconstructed MC samples, the same event and track selection criteria, as for real data, were used. The response of the detector was simulated based on the Geant package [31] (version 3.21), so the limited acceptance of the NA61/SHINE detector was also included in the reconstructed MC samples used to prepare the T MC res templates. Both the template and the data histograms were computed in selected bins of K * (892) 0 rapidity y (calculated in the center-of-mass reference system) and transverse momentum p T .
The T MC res and T DAT A mix histograms in the fit function given by Eq. (1) were normalized to have the same numbers of pairs as the real data histogram in the invariant mass range from 0.6 to 1.6 GeV. The symbols a, b and c in Eq. (1) are the normalization parameters of the fit (a + b + c = 1). They describe the contributions of T MC res , T DAT A mix and BW to the invariant mass spectra. The mass and width of the K * (892) 0 are the parameters of the Breit-Wigner shape obtained within the mass window m 0 ± 4Γ 0 . The values received from total fit 2 (see Fig. 3 or 4 (right)) were used to obtain the uncorrected numbers of K * (892) 0 mesons (the section below). For each m K + π − invariant mass bin in Fig. 3 and 4 (right), the bin content N bin (m K + π − ) was calculated as:  where N raw (m K + π − ) is the raw production in a given m K + π − bin, and a, b, T MC res (m K + π − ) and T DAT A mix (m K + π − ) are described in Eq. (1). The statistical uncertainty of N bin (m K + π − ) can be expressed as (the notation (m K + π − ) is omitted for simplifying the presentation of the formula): where ∆N raw , ∆T MC res and ∆T DAT A mix are the standard statistical uncertainties taken as the square root of the number of entries. For T MC res and T DAT A mix histograms the number of entries had to be properly normalized. Due to high statistics of Monte Carlo and mixed events, the uncertainties of parameters a and b were neglected.
In order to subtract a possible residual background (red curves) in Figs. 3 and 4 (right) (it looks negligible in these (y,p T ) intervals but is more significant in others), a fit of the blue histograms was performed as the last step using the function given by Eq. (5): where d, e, f , and g are free parameters of the fit, and the Breit-Wigner (BW) component was described by Eq. (2). The results are presented in Figs. 3 and 4 (right). The red lines here (polynomial background) illustrate the remaining residual background (Eq. (5) without BW component) and the brown curves (total fit 2) the sum of residual background and BW signal distribution (Eq. (5)). In the end, the uncorrected number of K * (892) 0 resonances (for each separate rapidity and transverse momentum bin) was obtained as the integral (divided by the bin width) over the BW signal of total fit 2 in Figs. 3 and 4 (right). The integral was calculated in the mass window m 0 ± 4Γ 0 .
3.7. Uncorrected numbers of K * (892) 0 Table 2 presents the uncorrected numbers of K * (892) 0 mesons, N K * (y, p T ), as obtained from the extraction procedure described in Sec. 3.6. The values are shown with statistical uncertainties. Due to limited statistics of data two kinds of binning were proposed: (i) one large bin/range of transverse momentum (0 < p T < 1.5 GeV/c) and two (40 GeV/c) or four (80 GeV/c) bins in rapidity (upper part of Table 2), (ii) one large bin/range of rapidity (0 < y < 1.5) and four bins in transverse momentum (lower part of Table 2). Binning presented in the upper part of Table 2 was used to obtain rapidity spectra of K * (892) 0 mesons, whereas binning illustrated in the lower part of Table 2 was used to compute transverse momentum and transverse mass spectra, as well as the p T dependence of the fitted mass and width of the K * (892) 0 resonance.
For each bin of (y, p T ) in Table 2 Table 2.: The uncorrected numbers of K * (892) 0 mesons, N K * (y, p T ), obtained from the extraction procedure described in Sec. 3.6 for inelastic p+p interactions at 40 GeV/c (middle column) and 80 GeV/c (right column). The values are shown with statistical uncertainties. Upper part: binning used to obtain y spectra of K * (892) 0 (see Fig. 8).
Lower part: binning used to obtain p T and m T spectra of K * (892) 0 , as well as the p T dependence of m K * and Γ K * (see Figs. 6, 7, 5).

Correction factors
The procedure of determining the uncorrected numbers of K * (892) 0 mesons was described in Sec. 3.7. These numbers need to be corrected for the effects such as identification inefficiency, geometrical acceptance, track and event reconstruction inefficiencies, and losses of inelastic p+p events due to the trigger bias (S4). In order to obtain the corrected numbers of K * (892) 0 mesons, produced in inelastic p+p interactions, two corrections were applied to the extracted raw numbers of K * (892) 0 resonances: (i) The loss of the K * (892) 0 mesons due to the dE/dx requirements was corrected by a constant factor: where K + = 0.84900 and π − = 0.99605 are the probabilities (based on the cumulative Gaussian distribution) for K + or π − to lie within (−1.2σ; +1.8σ) or (−2.7σ; +3.3σ) around the empirical parametrization of Bethe-Bloch value.
(ii) The losses due to geometrical acceptance, reconstruction efficiency, trigger bias (S4), detector acceptance as well as the quality cuts applied in the analysis were corrected with the help of a detailed Monte Carlo simulation. In the MC samples, the width of the K * (892) 0 resonance was simulated according to the known PDG value [33]. The correction factors were based on 19.7 × 10 6 (p beam = 40 GeV/c) and 19.8 × 10 6 (p beam = 80 GeV/c) inelastic p+p events produced by the Epos1.99 event generator [29]. The validity of these events for calculation of the corrections was verified in Refs. [23,34]. The particles produced in the generated events were tracked through the NA61/SHINE apparatus using the Geant package [31] (version 3.21). As the next step, the TPC response was simulated by dedicated NA61/SHINE software packages, which take into account all known detector effects. Then, the simulated events were reconstructed with the same software as used for the real data. Finally, the same selection cuts were applied (with the exception of the identification cut: dE / dx versus total momentum p lab ; instead, the matching procedure between reconstructed and simulated tracks was applied -see below).
For a given y and p T bin, the correction factor c MC (y, p T ) was computed as: where: -N gen K * (y, p T ) denotes the number of K * (892) 0 mesons (that decay into K + π − pairs) generated in a given (y,p T ) bin, -N sel K * (y, p T ) denotes the number of K * (892) 0 mesons (that decay into K + π − pairs) reconstructed and selected by the cuts in a given (y, p T ) bin. In this analysis the reconstructed charged particles were matched to the simulated K + and π − mesons based on the number of clusters and their positions. Then the invariant mass was calculated for all K + π − pairs. The reconstructed number of K * (892) 0 resonances was obtained by repeating the same steps (template method) as in raw experimental data (details are described in Sec The statistical uncertainty of c MC (y, p T ) was calculated assuming that N sel K * (y, p T ) is a subset of N gen K * (y, p T ) and the uncertainty of their ratio is governed by a binomial distribution. The uncertainty originating from the N sel events /N gen events ratio was found to be negligible. The final uncertainty of c MC (y, p T ) was then calculated as follows: The obtained values of correction factors c MC (y, p T ), together with statistical uncertainties, are presented in Table 3 for all considered (y, p T ) bins. 1.079 ± 0.022 2.86 ± 0.12 Table 3.: The correction factors c MC (y, p T ) with statistical uncertainties for 40 GeV/c (middle column) and 80 GeV/c (right column). Upper part: binning used to obtain y spectra of K * (892) 0 (see Fig. 8). Lower part: binning used to obtain p T and m T spectra of K * (892) 0 , as well as the p T dependence of m K * and Γ K * (see Figs. 6, 7, 5).

Statistical and systematic uncertainties
The statistical uncertainties of the corrected double-differential K * (892) 0 yields (see Eq. (9)) include the statistical uncertainties of the correction factor c MC (y, p T ) (see Eq. (8)) and the statistical uncertainties ∆N K * (y, p T ) (see Sec. 3.7) of the uncorrected number of K * (892) 0 resonances. The correction c dE / dx has no statistical uncertainty. The final expression for statistical uncertainty reads: The uncorrected numbers of K * (892) 0 mesons (and later on the corrected yields), the K * (892) 0 mass and width parameters, and other quantities depend on the details of signal extraction procedure and the event and track quality cuts. These two groups of effects were studied in order to estimate the systematic uncertainties.
(I) The uncertainties estimated by changing the signal extraction procedure: (iii) the dE / dx cuts, (−1.2σ; +1.8σ) for K + and (−2.7σ; +3.3σ) for π − , were changed into (−0.7σ; +1.3σ) for K + and (−2.2σ; +2.8σ) for π − (narrower cut), as well as (−1.7σ; +2.3σ) for K + and (−3.2σ; +3.8σ) for π − (wider cut), (iv) the minimum required total number of points in all TPCs for K + and π − candidates was modified from 30 to 25 and 35, (v) the minimum required number of clusters in both VTPCs for K + and π − candidates was modified from 15 to 12 and 18, (vi) the impact parameter (distance between the extrapolated track and the interaction point) cuts for the tracks were turned off.
(III) The uncertainties due to the limited precision of magnetic field calibration.
The NA61/SHINE magnetic field strength was verified with a precision of better than 1% by studying the K 0 S and Λ invariant mass distributions [36]. As in the previous paper [9] in order to test how the magnetic field calibration influences the results, the momentum components of K * (892) 0 decay products (K + and π − ) were varied by ±1%.
For each of the possible sources described above, the partial systematic uncertainty ∆ sys,i was conservatively determined as half of the difference between the lowest and the highest value obtained by varying the given parameter (statistical uncertainties were not considered while evaluating systematic uncertain-   Table 4. Within uncertainties, the values of Γ K * for both studied beam momenta (40 GeV/c and 80 GeV/c) are consistent with the PDG reference value (dashed horizontal line in Fig. 5 (bottom)). For 40 GeV/c data, the m K * values are also in agreement with the PDG reference value (dashed horizontal line in Fig. 5 (top)). For 80 GeV/c beam momentum, the observed m K * values seem to be slightly smaller than the reference value provided by the PDG. The comparisons of NA61/SHINE mass and width of K * (892) 0 resonances with STAR p+p results are shown in Sec. 5. The numerical data are listed in Table 4. The dashed horizontal lines represent PDG values m 0 = 895.55 MeV and Γ 0 = 47.3 MeV [28]. For a comparison the previous NA61/SHINE results [9] for p+p interactions at 158 GeV/c ( √ s NN = 17.3 GeV) are also shown (they were obtained in 0 < y < 0.5).

Double-differential K * (892) 0 spectra
The double-differential yields d 2 n dy d p T of K * (892) 0 mesons in inelastic p+p interactions at 40 GeV/c and 80 GeV/c were computed from Eq. (9). They are presented in Fig. 6 in bins of transverse momentum (see Sec. 4.3). The d 2 n dy d p T values in bins of rapidity were used to obtain the dn dy spectra presented in Fig. 8 Table 4.: The numerical values of mass and width of K * (892) 0 mesons fitted in 0 < y < 1.5 and presented in Fig. 5. The first uncertainty is statistical, while the second one is systematic.
4.3. K * (892) 0 transverse momentum and transverse mass spectra Figure 6 presents the double-differential yields of K * (892) 0 mesons as function of p T for rapidity range 0 < y < 1.5. The corresponding numerical values are listed in Table 5. In order to determine the inverse slope parameter T of transverse momentum spectra the function: was fitted to the measured data points shown in Fig. 6. The parameter A represents the normalization factor. The inverse slope parameters, resulting from the fits, are quoted in the figure legends.   Table 5.: The numerical values of double-differential yields d 2 n dy d p T presented in Fig. 6, given in units of 10 −3 (GeV/c) −1 . The first uncertainty is statistical, while the second one is systematic.
The transverse mass (m T ≡ p 2 T + m 2 0 ) spectra 1 m T d 2 n dm T dy were obtained based on d 2 n dy d p T spectra according to the relation: The results are presented in Fig. 7, together with the previous NA61/SHINE measurement at 158 GeV/c [9].
The numerical values for this analysis are displayed in Table 6. At higher energies the m T spectra seem to exhibit the concave shape with respect to the fitted exponential parametrization.   Table 6. The solid lines represent function given by Eqs. (11) and (12) with A and T parameters taken from Fig. 6. For a comparison the previous NA61/SHINE results [9] for p+p interactions at 158 GeV/c are also shown (they were obtained in 0 < y < 0.5).  Table 6.: The numerical values of double-differential yields 1 m T d 2 n dm T dy given in units of 10 −3 (GeV) −2 and presented in Fig. 7 (for 40 GeV/c and 80 GeV/c data); the values of m T − m 0 specify the bin centers. The first uncertainty is statistical, while the second one is systematic.
The inverse slope parameters of transverse momentum spectra (Fig. 6) in 0 < y < 1.5 were found to be T = (153 ± 29 ± 13) MeV for p beam = 40 GeV/c and T = (153 ± 30 ± 9) MeV for p beam = 80 GeV/c. The statistical uncertainty (the first one) is equal to the uncertainty of the fit parameter, and the systematic uncertainty was estimated in the way described in Sec. 3.10. In the previous analysis of NA61/SHINE the value of T = (173 ± 3 ± 9) MeV was obtained in 0 < y < 0.5 for p+p interactions at p beam = 158 GeV/c [9] (see Fig. 7). Finally, also in p+p collisions at 158 GeV/c the NA49 experiment measured the T parameter of the p T spectrum (for rapidity range 0.2 < y < 0.7) and published a value T = (166 ± 11 ± 10) MeV [10].

K * (892) 0 rapidity spectra
The K * (892) 0 rapidity distributions dn dy , presented in this paper, were obtained in transverse momentum range 0 < p T < 1.5 GeV/c. They were computed from d 2 n dy d p T values (in rapidity bins) multiplied by the width of the transverse momentum bin (1.5). The uncertainties were also obtained by multiplying the uncertainties of d 2 n dy d p T by 1.5. The spectra are presented in Fig. 8 together with the previous NA61/SHINE results obtained in the full p T range for p+p interactions at 158 GeV/c [9]. The numerical values for this analysis are displayed in Table 7. The data points presented in Fig. 8 were fitted with a Gaussian function: that allowed to determine the width σ y of the K * (892) 0 rapidity distribution. The parameter A represents the normalization factor. Note that in the fit function the mean value of the Gaussian shape was fixed at y = 0. The fit parameters were also used to compute the mean multiplicity K * (892) 0 (details are given in Sec. 4.5). The statistical uncertainty of σ y was taken from the fit, and the systematic one was estimated in the way described in Sec. 3.10. The numerical values of σ y and K * (892) 0 are shown in Table 7.   Table 7. The solid lines represent the function given by Eq. (13). For comparison the previous NA61/SHINE results [9] for p+p interactions at 158 GeV/c are also shown (they were obtained in the full transverse momentum range; p T -integrated and extrapolated rapidity spectrum). For 158 GeV/c the first (light blue) point (y < 0) was not included in the fit (see Ref. [9] for details).

Mean multiplicity of K * (892) 0 mesons
The mean multiplicities of K * (892) 0 mesons were obtained based on rapidity distributions presented in Fig. 8. Assuming rapidity symmetry around y = 0, the mean multiplicity K * (892) 0 was calculated as the sum of measured points in Fig. 8 Table 7.: The numerical values of rapidity distributions presented in Fig. 8. The first uncertainty is statistical, while the second one is systematic. Additionally, the table shows the widths of the Gaussian fits to the dn dy distributions and the mean multiplicities of K * (892) 0 mesons (see Sec. 4.5 of the text for details). unmeasured region: where for 80 GeV/c data: and for 40 GeV/c data: The function f (y) is described by Eq. (13). The statistical uncertainty of K * (892) 0 was determined as: where ∆ dn dy is the statistical uncertainty of dn dy point and ∆y is the rapidity bin width (equal 0.5 for each of the i-th dn dy points). The systematic uncertainty of K * (892) 0 was estimated in the way described in Sec. 3.10. The results are listed in Table 7 and presented in Fig. 9. The mean multiplicities of K * (892) 0 mesons in inelastic p+p collisions were found to be (35.1 ± 1.3(stat) ± 3.6(sys)) · 10 −3 at 40 GeV/c and (58.3 ± 1.9(stat) ± 4.9(sys)) · 10 −3 at 80 GeV/c.

Comparison with world data and model predictions
Comparisons of the NA61/SHINE measurements with publicly available world data are presented. The results are also confronted with predictions of Epos1.99 and statistical models.  [9] for p+p interactions at 158 GeV/c is also shown. The mean multiplicities were obtained for 0 < p T < 1.5 GeV/c (two lower energies, this analysis) or for the full phase space [9] (the highest energy). The vertical bars represent the total uncertainties (square root of the sum of squares of statistical and systematic uncertainties). The dashed line is added as a guide to the eye.

Mass and width of K * (892) 0
In Fig. 10 the results of NA61/SHINE for K * (892) 0 mass and width in inelastic p+p collisions (this analysis and Ref. [9]) are compared to p+p results from STAR at RHIC and the PDG reference values (for STAR the mass and width of K * 0 meson peak were calculated as the averaged measurements from K * (892) 0 and K * (892) 0 invariant mass spectra). Similar plots presenting Pb+Pb and Au+Au results (NA49, ALICE, STAR) can be found in Ref. [9].
The obtained NA61/SHINE measurements of m K * and Γ K * are close to the PDG reference values. However, somehow lower m K * values may be seen for 80 GeV/c p+p data. For p+p collisions at RHIC energy, the STAR experiment also measured lower K * 0 mass, especially at lower transverse momenta.

Comparison of results with Epos1.99 predictions and NA49 measurements
The NA61/SHINE results on rapidity spectra and mean multiplicities were compared to predictions of the Epos1.99 [29] model of hadron production. The rapidity spectra are presented in Fig. 11, and the numerical values of mean multiplicities are listed in Table 8. For comparison, the previous NA61/SHINE result for 158 GeV/c [9] is also included in the table (it was obtained from p T -intergated and extrapolated dn dy spectrum, thus resulting in K * (892) 0 measured in the full phase space [9]). It can be seen that the Epos1.99 model overestimates K * (892) 0 production in inelastic p+p collisions at all three SPS beam momenta.   Figure 10.: The transverse momentum dependence of mass and width of K * (892) 0 (or K * 0 ) mesons obtained in p+p collisions by NA61/SHINE (this analysis and Ref. [9]) and STAR [11]. For STAR the averaged (K * 0 ) measurements of K * (892) 0 and K * (892) 0 are presented. The horizontal lines represent PDG values [28]. Table 8 also includes the comparison of p+p results for 158 GeV/c with NA49 [10]. The NA49 experiment used one wide p T bin (0 < p T < 1.5 GeV/c; similarly to the NA61/SHINE analysis of 40 and 80 GeV/c data) and the K * (892) 0 was obtained from the dn dy spectrum as the integral under the Gaussian function in the range −3 < y < 3 [35]. Within the estimated uncertainties, the results of both experiments were consistent.
paper, the measured NA61/SHINE K * (892) 0 multiplicities are compared with predictions [37] of the HRG model with parameters obtained by fitting the NA61/SHINE p+p data. Figure 12 presents the energy dependence of K * (892) 0 to K * (892) 0 HRG ratios for the HRG model [37] in the Canonical Ensemble (CE). The upside-down red triangles correspond to the HRG fits with the φ meson multiplicities included, whereas violet triangles represent the situation where the φ mesons were not included in the HRG model fits. Additionally, the NA61/SHINE p+p point at 158 GeV/c was compared to the HRG model prediction within the Grand Canonical Ensemble (GCE) formulation [37,38] (blue star symbol in Fig. 12). In Fig. 12 the total uncertainty of K * (892) 0 was taken as the square root of the sum of squares of statistical and systematic uncertainties. The uncertainty of the K * (892) 0 to K * (892) 0 HRG ratio (vertical axis) was taken as the final uncertainty of K * (892) 0 divided by K * (892) 0 HRG . The Hadron Resonance Gas model in the CE agrees with the NA61/SHINE p+p data at p beam = 40-158 GeV/c but only when the φ meson is excluded from the fit. The Authors of Ref. [37] stress that the inclusion of the φ meson multiplicities in thermal fits significantly worsens the HRG model fit quality. But surprisingly, the GCE statistical model well describes the K * (892) 0 yield in the small p+p system (point for 158 GeV/c). Note that the K * /K ratios in p+p collisions at higher energies are also consistent with the GCE statistical model predictions [14,15,19,39]. The numerical values used to prepare Fig. 12 are presented in Table 9 of Appendix A.  [37] for the fit with φ mesons included (upside-down triangles) and the fit with φ meson excluded (triangles). The star shows the K * (892) 0 [9] divided by the HRG model prediction for the Grand Canonical Ensemble [37,38]. The numerical values of K * (892) 0 and K * (892) 0 HRG are listed in Table 9 of Appendix A.

K * (892) 0 over charged kaon ratios
The system size dependence or multiplicity dependence of K * to charged kaon ratios may allow estimating the time interval between chemical and kinetic freeze-out in nucleus-nucleus (A+A) collisions [8]. This is done based on the ratio of the K * /K produced in A+A and p+p collisions. The K * (892) 0 / K + and K * (892) 0 / K − ratios in p+p are shown in Fig. 13, and the corresponding numerical values are listed in Table 10 of Appendix A. Together with future NA61/SHINE measurements in Be+Be, Ar+Sc, and Xe+La collisions, it will allow estimating the time between freeze-outs for these nucleus-nucleus systems at three SPS energies.  Table 10 of Appendix A (p+p at 40 and 80 GeV/c) and in Ref. [9] (NA61/SHINE p+p data at 158 GeV/c).

Blast-Wave model fits
The fits within the Blast-Wave models allow obtaining thermal freeze-out temperature (T f o ) and transverse flow velocity (β T ) of the system. The transverse mass spectra of K * (892) 0 mesons (this analysis and Ref. [9]) and other particles previously reported by NA61/SHINE (charged pions, charged kaons, protons, anti-protons [40], φ mesons [41]) were fitted within the Blast-Wave model [6] with β T independent of the radial position in the thermal source. The fitted formula follows: where I 0 and K 1 are the modified Bessel functions, A i are the fitted normalization parameters, and index i refers to different particle species. The fit parameter ρ is related to the transverse flow velocity by ρ = tanh −1 β T . The results of a simultaneous fit to the m T distributions of different particle species are presented in Fig. 14 [42][43][44] at the same beam momenta.  Figure 14.: The transverse mass spectra of K * (892) 0 mesons (0 < y < 1.5 for 40 and 80 GeV/c from this analysis, and 0 < y < 0.5 for 158 GeV/c from Ref. [9]) and other hadrons previously measured by NA61/SHINE (charged pions, charged kaons, protons, anti-protons [40] in 0 < y < 0.2, and φ mesons [41] in 0 < y < 0.3; for π − at 80 GeV/c the rapidity range 0.2 < y < 0.4 was used instead of 0 < y < 0.2) fitted within the BW model [6] described by Eq. (18). Results for 40 GeV/c (left), 80 GeV/c (middle), and 158 GeV/c (right) beam momenta. For all points the vertical uncertainty bars represent total uncertainties (square root of the sum of squares of statistical and systematic uncertainties). The fits were performed in the range 0 < m T − m 0 < 1 GeV. The resulting fit parameters are displayed in the legends.

Summary
This publication presents the NA61/SHINE measurements of K * (892) 0 meson production via its K + π − decay mode. The results were obtained for inelastic p+p collisions at beam momenta 40 GeV/c and 80 GeV/c ( √ s NN = 8.8 and 12.3 GeV). The template method was used to extract raw K * (892) 0 signals. In this method, the background is described as a sum of two contributions: background due to uncorrelated pairs modeled by event mixing and background of correlated pairs modeled by Epos1.99.
The fits to background-subtracted invariant mass spectra were used to obtain the masses and widths of the K * (892) 0 resonance. The NA61/SHINE values, for different transverse momentum bins, are generally close to the PDG results, however, a small deviation from the reference value may be observed for K * (892) 0 mass at 80 GeV/c.
The NA61/SHINE results were compared with predictions of the Epos1.99 model and the Hadron Resonance Gas model. Epos1.99 overestimates K * (892) 0 production in p+p collisions at the SPS energies. The Canonical Ensemble formulation of the HRG model gives a good description of p+p data provided that the φ meson is excluded from the fits.