Measurement of the central exclusive production of charged particle pairs in proton-proton collisions at $\sqrt{s} = 200$ GeV with the STAR detector at RHIC

We report on the measurement of the Central Exclusive Production of charged particle pairs $h^{+}h^{-}$ ($h = \pi, K, p$) with the STAR detector at RHIC in proton-proton collisions at $\sqrt{s} = 200$ GeV. The charged particle pairs produced in the reaction $pp\to p^\prime+h^{+}h^{-}+p^\prime$ are reconstructed from the tracks in the central detector, while the forward-scattered protons are measured in the Roman Pot system. Differential cross sections are measured in the fiducial region, which roughly corresponds to the square of the four-momentum transfers at the proton vertices in the range $0.04~\mbox{GeV}^2<-t_1 , -t_2<0.2~\mbox{GeV}^2$, invariant masses of the charged particle pairs up to a few GeV and pseudorapidities of the centrally-produced hadrons in the range $|\eta|<0.7$. The measured cross sections are compared to phenomenological predictions based on the Double Pomeron Exchange (DPE) model. Structures observed in the mass spectra of $\pi^{+}\pi^{-}$ and $K^{+}K^{-}$ pairs are consistent with the DPE model, while angular distributions of pions suggest a dominant spin-0 contribution to $\pi^{+}\pi^{-}$ production. The fiducial $\pi^+\pi^-$ cross section is extrapolated to the Lorentz-invariant region, which allows decomposition of the invariant mass spectrum into continuum and resonant contributions. The extrapolated cross section is well described by the continuum production and at least three resonances, the $f_0(980)$, $f_2(1270)$ and $f_0(1500)$, with a possible small contribution from the $f_0(1370)$. Fits to the extrapolated differential cross section as a function of $t_1$ and $t_2$ enable extraction of the exponential slope parameters in several bins of the invariant mass of $\pi^+\pi^-$ pairs. These parameters are sensitive to the size of the interaction region.


Introduction
The study of exclusive production of meson and baryon pairs has long been recognised as an important ground for Quantum Chromodynamics (QCD) tests. Exclusive production of pion and kaon pairs has been studied both theoretically [1][2][3] and experimentally, in two-photon collisions at lepton colliders [4][5][6][7] and via photoproduction [8][9][10][11] and deep inelastic scattering [12,13] in lepton-proton and heavy-ion experiments. Exclusive production of meson and baryon pairs belongs to the class of Central Exclusive Production (CEP) processes. In hadron-hadron collisions, CEP processes can proceed via Double Pomeron Exchange (DPE), photon-Pomeron exchange or photon-photon exchange, where the Pomeron is a colour-singlet object with internal quantum numbers of the vacuum, see, e.g., [14,15]. Although several properties of diffractive scattering at high energies are described by the phenomenology of Pomeron exchange in the context of Regge theory [16], the exact nature of the Pomeron still remains elusive. This paper presents a measurement of CEP of π + π − , K + K − and pp pairs in pp collisions at a centre-of-mass energy of √ s = 200 GeV with the STAR detector at RHIC. Differential and integrated cross sections are measured in a fiducial region and compared to phenomenological predictions. The fiducial region roughly corresponds to the square of the fourmomentum transfers at the proton vertices in the range 0.04 GeV 2 < −t 1 , −t 2 < 0.2 GeV 2 and invariant masses of the charged particle pairs up to a few GeV. Throughout the paper the convention c = = 1 is used.

Theoretical framework and current experimental situation
Over the last decade, one could observe a renewal of interest in studies of CEP processes in high energy proton-(anti)proton collisions (see, e.g., [17][18][19] for review and further references). CEP processes in hadron-hadron collisions provide an especially clean environment to study the nature and quantum numbers (spin, parity) of centrally-produced resonance states [18]. In proton-proton collisions, the CEP reaction may be written in the form pp → p ⊕ X ⊕ p , (2.1) where the ⊕ symbols denote the presence of large rapidity gaps which separate the final state system X from the diffractively-scattered protons. This process in the DPE mode and with the hadronic final state, X, consisting of just an oppositely-charged particle pair is schematically shown in Fig. 1(left). In Fig. 1(right), its representation within perturbative QCD (pQCD) is shown in the two-gluon approximation. The scattered protons emerge intact from the collision at small polar angles with respect to the incoming beams and can be detected with special tagging devices. The final state, X, can be fully measured at central rapidities. The upper limit of the invariant mass, M X , of the system X depends on the energy of the colliding particles. At the LHC, this upper limit can reach above 100 GeV and the central production of even the Higgs or supersymmetric particles might be possible [20]. Recently, CEP of pion pairs in pp collisions at √ s = 5.02, 7 and 13 TeV has been reported by the CMS experiment [21,22]. The LHCb experiment has measured photo-induced CEP of J/ψ and ψ(2S) states in pp collisions at √ s = 13 TeV [23]. At the Tevatron, it was possible to study CEP of π + π − pairs with invariant masses up to a few GeV in pp collisions at several centre-of-mass energies up to √ s = 1.96 TeV [17,24]. The above measurements were performed without forward-proton tagging. The experiments at CERN at the ISR [25] and the SPS [26] have provided measurements of many CEP processes with forward proton tagging, however at significantly smaller centre-of-mass energies (62 GeV for ISR and 30 GeV for SPS).
In the present analysis, X stands for either continuum or resonant production of π + π − , K + K − or pp pairs in the non-perturbative regime up to around M X = 3 GeV. At the relatively high centre-of-mass energy of √ s = 200 GeV, only small contributions from Reggeon exchange are expected. Contributions from photon-Pomeron and photon-photon processes are also not significant and are additionally suppressed at −t > 0.04 GeV 2 . Therefore, the DPE process is expected to be dominant. The DPE process can be regarded as a spin and parity filter, i.e., the system X must have even spin and positive parity, so only exclusive production of f 0 and f 2 resonances are allowed on top of the continuum. In general, the resonance and continuum contributions may interfere. Calculations of the hadron mass spectrum in this domain were first done for only the continuum production [18,27,28] of π + π − or K + K − pairs using an approach based on Regge theory. In these models, the parameters of the Pomeron and sub-leading Reggeon exchanges were adjusted to describe the total and elastic πp or Kp scatterings. In this approach, the amplitude for the p + p → p + ππ(KK) + p process is expressed in terms of the product of two amplitudes describing the interaction of each of the two protons with one of the two mesons. The intermediate meson form factor is parameterised with one of three functions: an exponential, exp (t − m 2 h )/Λ 2 off , an Orear-like function, exp −b( √ −t + m 2 h + a 2 − a) , or a power-like function, 1/(1 − (t − m 2 h )/a 0 ), wheret is the square of the four momentum transfer at the Pomeron-meson vertex and Λ 2 off , b, a and a 0 are free parameters. The models can be supplemented to include absorption effects (shown symbolically on Fig. 1) which are needed to calculate the 'survival probability' for no additional soft re-scatterings between the colliding protons or the final state mesons. Absorption corrections are related to the non-perturbative interaction in the initial or final state of the reaction. The re-scattering leads to suppression of the cross section and distortion of the distributions of the kinematic variables. The suppression factor depends on the collision energy. They usually reduce the cross sections, even by a factor of 5 at RHIC energy and a factor of 10 at LHC energies [29]. Recently, the production of a variety of resonances: the f 0 (500), f 0 (980) and f 2 (1270) decaying to π + π − , the f 0 (980), f 0 (1500), f 0 (1710), f 2 (1270) and f 2 (1525) decaying to K + K − , and the f 0 (2020), f 0 (2100) and f 0 (2200) decaying to pp, were studied theoretically [30][31][32] including interference effects between resonant and non-resonant amplitudes. The calculations are based on a tensor Pomeron model. The amplitudes for the processes are

Experimental setup
The data used in this analysis were collected by the STAR experiment at RHIC [41] in 2015 in proton-proton collisions at √ s = 200 GeV and correspond to an integrated luminosity of 14.2 pb −1 . A detailed description of the STAR detector is given in Ref. [42].
The forward-scattered protons are measured in the Roman Pot (RP) system adopted from the pp2pp experiment [43] at RHIC. This is schematically shown in Fig. 2 below the beam-line, see Fig. 2 (bottom). The RPs are situated downstream of the DX dipole magnets responsible for head-on targeting of the incoming beams and for bending outgoing beams back into the respective accelerator pipelines. The constant and uniform magnetic field of the DX magnet works as a spectrometer, and thus knowledge of the scattered proton's trajectory allows reconstruction of its momentum. Each RP houses a package of 4 silicon strip detector planes -two with vertical and two with horizontal orientation of the strips -allowing measurement of the position of a proton in the transverse plane. The strip pitch is about 100 µm, resulting in a spatial resolution of about 30 µm. Scintillator counters placed inside each RP station allow for triggering on forward protons and also provide timing information with 0.5 ns resolution. Measurement of a pair of charged particles produced in the final state at central rapidity is performed using the Time Projection Chamber (TPC) [44], which provides tracking for charged particles in the 0.5 Tesla solenoidal magnetic field. The TPC covers the pseudorapidity range of |η| < 1.8 in full azimuthal angle. The TPC is used to determine the momenta of the charged particles, and also helps in locating the position of the collision vertex. The tracking efficiency is ∼ 85% for |η| < 1, but falls to 50% at |η| ∼ 1.3. The measurement of the specific energy loss in the TPC gas, dE/dx, is used for particle identification. Furthermore, to extend the particle identification power of the STAR detector performed by the TPC, a Time-Of-Flight (TOF) detector [45] is placed around the TPC covering a pseudorapidity range of |η| < 0.9. The TOF detector is a system of adjacent Multi-gap Resistive Plate Chambers. In addition to being a precise timing detector, it is used to measure event multiplicity at the trigger level and to discriminate TPC tracks arrived in preceding/posterior bunch crossings (out-of-time pile-up) from the in-time tracks. Using the TOF timing information together with the momentum and path length reconstructed in the TPC allows for the particle mass determination.
For the TOF efficiency study, an unbiased sample with tracks reconstructed using Heavy Flavor Tracker (HFT) [46] hits in addition to TPC hits was analysed. This sample provides a clean source of in-time tracks. The HFT is a system of multi-layer silicon pixel and strip detectors. It improves the impact parameter resolution of the STAR tracking system and enables reconstruction of secondary decay vertices of open heavy flavour hadrons.
To suppress non-exclusive background, the Beam-Beam Counters (BBCs) [47] and the Zero-Degree Calorimeters (ZDCs) [48] are used. The BBCs are scintillator detectors placed in the endcap regions of the STAR detector and covering pseudorapidity ranges of 2.1 < |η| < 3.3 (large BBC tiles) and 3.3 < |η| < 5 (small BBC tiles). The ZDCs are used to tag neutral particles which leave the interaction region close to the beam direction.

Event reconstruction
In order to reconstruct the position and momentum of the scattered protons, a clustering procedure is applied for each RP detector plane separately. A cluster is formed by a continuous series of Si strips with signals well above the pedestal. Pairs of matched clusters found in detector planes measuring the same coordinate define the (x, y) coordinates of space points for a given RP. Correlating space points between RP stations, and reconstructing the proton kinematics, relies on an alignment procedure, which is carried out using elastic scattering data for each run separately, as described in Ref. [49]. A track is formed based on one or two points reconstructed in the two RP detector stations on the same side of the IP. Using elastic scattering events reconstructed with tracks formed from the two points, the average transverse position of the primary vertex was measured to be x IP = 0.42±0.04 mm and y IP = 0.46 ± 0.05 mm. With the average transverse position of the vertex and twopoint proton tracks, the proton transverse momentum (p x , p y ) can be reconstructed, and hence the value of the Mandelstam variable t. For single point proton tracks, the transverse momentum is reconstructed assuming that the scattered proton energy is equal to the beam energy. Such an approximation is justified because the proton loses on average less than 1% of its initial energy for events with M X < 2 GeV.
Particle pair identification is performed using the combined information from the TPC and TOF detectors for both tracks simultaneously. The compatibility of the track's dE/dx with that expected for a given particle (h = π, K, p) is determined using the quantity where (dE/dx) h is the Bichsel [50] expectation for particle type h and σ is the relative resolution of dE/dx for a given track. From nσ h for each of the two tracks, the χ 2 statistic for a hh pair hypothesis is calculated: The time at which a particle is detected in the TOF system is used to reconstruct its squared mass m 2 TOF . For this purpose, the time of the primary interaction is required; however, it is not known for CEP events. Instead, the unknown time of the primary interaction can be eliminated by assuming that both tracks present in an event are of the same type. In that case, the measured TOF time difference between particles is given by ∆t = L 1 1 + m 2 TOF /p 2 1 − L 2 1 + m 2 TOF /p 2 2 , where p 1,2 are tracks' momenta and L 1,2 are the lengths of the helical paths between the primary vertex and TOF hit associated with them. m 2 TOF can be then calculated per event. If m 2 TOF is negative, due to the detector's resolution, then it is set to zero.

Monte Carlo simulation
Monte Carlo (MC) simulations are used for the modelling of background contributions, unfolding of detector effects, calculation of systematic uncertainties and comparisons of models with the hadron-level cross section measurements.
In the DiMe event generator, four models for absorption are available. The prediction from "model 1", which is most consistent with data, is used in this analysis. DiMe predictions are also sensitive to the choice of meson form factor. Three different parameterisations of meson form factor are implemented. We chose an exponential form using Λ 2 off = 1.0 GeV 2 , which fits the present data better. A larger value of Λ 2 off was used [18] to fit DiMe predictions to ISR data, however the 50% normalisation uncertainty of the ISR data does not exclude Λ 2 off = 1.0 GeV 2 . In GenEx, the absorption corrections are not taken into account. However, the model developers estimated the suppression factor to be of the order of 2 − 5 (π + π − ) and 2 (K + K − ) [29]. To account for absorption, the π + π − cross sections obtained from GenEx are scaled by 0.25 to fit DiMe predictions for masses above 0.8 GeV, while the K + K − cross sections from GenEx are scaled by 0.45 to fit DiMe predictions for masses above 1.2 GeV. Above these limits, the absorption effects only weakly depend on pair mass. In the GenEx generator, we also use an exponential form for the meson form factor using Λ 2 off = 1.0 GeV 2 . Therefore the differences between GenEx and DiMe are almost entirely due to the absorption effects.
The MBR model [52] implemented in PYTHIA8 [53] was tuned to describe the inclusive cross section for central diffraction (CD), p + p → p + X + p , measured by the CDF experiment. In this model, the exclusive h + h − final state occurs from fragmentation and hadronisation of the central state based on the Lund string model. The MBR model implemented in PYTHIA8.244 allows generation of the central state starting from the mass threshold of 0.5 GeV, however the suggested value is 1.2 GeV. Therefore, PYTHIA8 expectations for very low masses are in question, but are shown for completeness. The obtained cross sections from PYTHIA8 are scaled by an arbitrary value of 0.25 for easier comparison with the data.
Single particle MC embedded into zero-bias data events were used to calculate the TPC and TOF reconstruction and matching efficiencies separately for: π + , π − , K + , K − , p and p. The GenEx sample of p + p → p + π + π − + p embedded into zero-bias data events was used to calculate the RP reconstruction efficiency for forward-scattered protons and also for closure tests of the full analysis chain. The inclusive CD and Minimum Bias (in-elastic) PYTHIA8 samples without embedding were used for calculation of the impact of systematic uncertainties in the subtraction of non-exclusive background, as well as for the comparisons with the data. Prior to the embedding, MC samples were passed through a detailed GEANT3 [54] simulation of the STAR central detector and the GEANT4 [55] simulation of the beam optics and RP detectors. All MC samples were then subjected to the same reconstruction and analysis software as applied to the data.
A fast and simplified MC simulation was used for estimation of the central pair particle identification (PID) efficiency and misidentification probability. In this simulation, the dE/dx and the times of detection of particles in the TOF detector were generated according to parameterisations obtained from the inclusive data and the full TPC/TOF simulations, while the amount of exclusive π + π − , K + K − and pp were chosen to describe the data.

Data sample and event selection
The CEP events were triggered by requiring signals in at least one RP station on each side of the IP, and at least 2 hits in the TOF to ensure the presence of at least two in-time tracks in the TPC. In addition, a lack of activity in both the small BBC tiles and the ZDC detectors is required to ensure the double gap in pseudorapidity topology characteristic of CEP events. In order to reduce pile-up events, or events involving proton dissociation, a veto is imposed on events containing signals in both upper and lower RP stations on the same side of the IP. 560 million CEP event candidates were triggered in total, corresponding to 14.2 pb −1 of integrated luminosity. An average trigger prescale of 5 was used during the entire data-taking period due to limited data acquisition bandwidth.
In the offline analysis, the protons tagged in the RPs are further required to have transverse momenta (p x , p y ) in the fiducial region defined as (6.1) This fiducial region is chosen to achieve high geometrical acceptance and track reconstruction efficiency and also to minimise systematic uncertainties. Figure 3 (left) shows the combined distributions of the momenta, p y vs. p x , of the diffractively scattered protons in exclusive h + h − events reconstructed with the East and West RP stations. The kinematic region used in the measurement, defined in Eq. (6.1), is enclosed with the black line. In this analysis, the CEP events must consist of only one pair of oppositely charged mid-rapidity particles, besides the two forward-scattered protons. Therefore, only events with exactly two opposite-sign TPC tracks, each matched with hits in the TOF and originating from a common vertex, are selected. These tracks have at least 25 hits (out of a possible 45). All tracks are required to be within the fiducial region defined by p track T > 0.2 GeV and |η track | < 0.7. The z-position of the event vertex obtained from the TPC tracks is limited to |z vtx | < 80 cm. The above limits on η track and z vtx were chosen to ensure high geometrical acceptance in the entire fiducial phase space. In addition, it is required that the z-position of the vertex obtained from the time difference of the signals from the forward protons in the RPs agrees with the TPC vertex within 36 cm, corresponding to three-and-a-half standard deviations. To further suppress the residual backgrounds, a veto is imposed on signals in the large BBC tiles (2.1 < |η| < 3.3), as well as on events with more than one additional TOF hit not matching either of the two TPC tracks. This mainly removes higher-multiplicity events where some particles are either not reconstructed in the TPC or produced outside the TPC acceptance. PID involves a few steps. First, the pp hypothesis is checked: 2) If the conditions (6.2) are satisfied, the pair is assumed to be pp. If not, the pair is checked for compatibility with the K + K − hypothesis: 3) If the pair is not compatible with either the pp or the K + K − hypothesis, it is assumed to be a π + π − pair if For pairs identified as kaons or protons, a more restrictive cut on each track's p T is imposed: p T > 0.3 GeV for K ± and p T > 0.4 GeV for p(p). In addition, it is required that the lower-p T track in the pair has p T < 0.7 GeV (K + K − ) or p T < 1.1 GeV (pp). These additional cuts are intended to constrain the fiducial range of high track reconstruction efficiency (lower cut) and high pair identification efficiency (upper cut). The criteria used for PID given by Eqs. (6.2), (6.3) and (6.4), and the min(p + T , p − T ) cut discussed above, were chosen to suppress exclusive background below 1% for π + π − and pp and below 3% for K + K − . The distributions of χ 2 dE/dx and m 2 TOF for all studied particle species are shown in . Similarly, the m 2 TOF distributions are shown after the corresponding χ 2 dE/dx cut. Finally, the missing transverse momentum in the event, p miss T , obtained from the transverse momenta of the protons tagged in the RPs and the tracks of the centrally produced pair, is required to be less than 75 MeV to suppress the non-exclusive background. Figure 5 (left column) shows the p miss T distributions for all studied particle species together with the p miss T distributions for like-sign control sample. After all the above selection cuts, the approximate numbers of CEP event candidates are 85600 for π + π − pairs, 930 for K + K − pairs and 70 for pp pairs in the final state.
Uncorrected invariant mass distributions of the π + π − , K + K − and pp pairs after the final selection cuts are shown in Fig. 5 (right column). The same-sign control sample is at the level of a few percent of the final sample.

Kinematic variables and fiducial region
The measurements are done differentially in several kinematic variables, which include: • the cosine of polar angle (cos θ CS ) and the azimuthal angle (φ CS ) of positively charged central particle in the Collins-Soper frame [56] 1 with typical resolutions of 0.005−0.01 and 1 • − 2 • , respectively.
The differential cross sections are obtained in the fiducial region defined by the kinematical cuts imposed on the forward-scattered protons given in Eq. (6.1), and by the cuts on the final state charged particles' pseudorapidities: |η| < 0.7 and transverse momenta p T > 0.2 GeV (π + π − ), p T > 0.3 GeV (K + K − ) and p T > 0.4 GeV (pp). In addition, in the case of K + K − and pp pairs, the fiducial volume is restricted to the region with the lower p T in the pair below 0.7 GeV or 1.1 GeV, respectively.      Figure 4. Distributions of χ 2 dE/dx (left column) and m 2 TOF (right column) for exclusive π + π − (top), K + K − (middle) and pp (bottom) candidates after final event selection. The dashed red line and arrow indicate the value of the cut imposed on the plotted quantity to select exclusive pairs of a given particle species. Yellow, blue and green histograms correspond to the fast exclusive MC simulation while magenta shows the estimated amount of non-exclusive background in the data.    (left column) and invariant mass of the charged particle pairs produced in the final state (right column) for π + π − (top), K + K − (middle) and pp (bottom) pairs. Distributions for opposite-sign and same-sign particle pairs are shown as black and red symbols, respectively. The vertical error bars represent statistical uncertainties. The horizontal bars represent bin sizes. Solid magenta histograms correspond to the estimated non-exclusive background, determined differentially from the number of counts in the hatched range 0.16 GeV < p miss T < 0.24 GeV, and extrapolated to the signal region indicated with dashed red line and arrow.

Background estimation
Background in the analysis arises from non-exclusive processes leading to correlated signals in the RP and TOF/TPC ('single source') and from coincidences of a signal in the RP with an uncorrelated signal in the TOF/TPC ('pile-up'). Other sources of background are exclusive processes in which the particle pair was misidentified.
The 'single source' non-exclusive contribution is dominated by Central Diffraction, where Y is any number of particles produced (but not measured) in addition to the measured h + h − pair. In the 'pile-up' background, the signal in the central detector almost always arises from an inelastic pp collision while the RP signal occurs due to 'pile-up' from real forward-going protons from elastic scattering, central diffraction, showering in single, double or non-diffractive events or beam-induced sources. All the above sources of background are estimated using a data-driven method. Both undetected particles in 'single source' events and the random character of 'pile-up' events lead to breaking the correlation between the central h + h − pair and the forward protons and to a much flatter p miss T distribution. This can be seen in Fig. 5 (left column), where the p miss T distribution starts to increase above 100 MeV. Background estimation is based on the extrapolation into the signal region of the second-degree polynomial function fitted to the signal-free region, as shown by the magenta histograms. The polynomial is constrained in the fit to vanish at p miss T = 0. This procedure is repeated differentially for all the kinematic variables presented later. As an example, the resulting background estimation is presented in Fig. 5 (right column) differentially in m(h + h − ) by the magenta histograms. The shape of the background as a function of p miss T is confirmed by MC predictions and like-sign events (shown by red points), which are, by definition, background. On average, this background amounts to 5.3% (π + π − ), 5.4% (K + K − ) and 12% (pp).
The exclusive background was estimated based on the fast MC simulation. The resulting distributions of χ 2 dE/dx and m 2 TOF for all particle species are presented in Fig. 4.

Corrections
For all cross section calculations, bin sizes are chosen to correspond to about three times the detector resolution so that migrations between bins could be neglected. Particles passing through the detector material lose some energy. To minimise biases from this effect, a correction procedure is applied during track-momentum reconstruction for both data and MC simulation based on the expected material budget for the given track. In this procedure, all tracks are assumed to be pions, therefore the reconstructed momenta of the remaining particle species exhibit some bias. For tracks identified as kaons and protons, an additional energy loss correction was applied based on the single particle MC. The correction is up to 10 and 20 MeV for low-momentum kaons and protons, respectively.
Several corrections were implemented to account for the limited efficiency of the measurement. The RP and TOF trigger efficiencies were estimated from the unbiased data to be ∼ 100% and 98%, respectively. The RP trigger efficiency was evaluated as the probability that the trigger was set when a proton was reconstructed in the given station. For events with exactly two TPC tracks, each matched with hits in the TOF and originating from a common vertex, and at most one additional TOF hit, the TOF trigger efficiency was estimated as the fraction of events that passed the TOF trigger conditions. The average number of inelastic collisions per bunch crossing varies between 0.2 and 0.9 and leads to a sizeable probability that the exclusive signal overlaps with another process, which causes signal loss due to the trigger or offline selections. This probability is estimated based on zero-bias data for each run independently and parameterised as a function of instantaneous luminosity for each of the four possible combinations of RPs topology. An event is rejected if the overlapping process produces a signal in the BBC, ZDC, TOF, or the RP station not belonging to the studied combination, since our selection criteria include vetos on these detectors. The overall veto efficiency varies between 40 − 80%. The efficiency of the |z vtx | < 80 cm selection cut was estimated for each RHIC fill independently based on the estimated values of the mean and standard deviation, assuming a normal distribution in the data. A typical vertex-cut efficiency is 88%.
The above corrections were applied independent of all kinematic variables, affecting only the normalisation. Corrections described below were applied as functions of the relevant variables, affecting the shapes of all measured distributions.
Protons successfully reconstructed in a RP station may still produce secondaries in the dead material in the other RP branch, which cause the trigger to veto the event. The probability to pass the veto trigger was estimated using the embedded GenEx sample as a 3D function of proton momenta p x , p y and z vtx . The same sample was used to study the proton reconstruction efficiency. Secondaries produced in the dead material or from 'pile-up' processes introduce inefficiencies in the reconstruction procedure. This inefficiency was also calculated using a data-driven tag-and-probe method using elastic scattering data, where one of the protons serves as a tag to probe the reconstruction efficiency of the second proton. The joint efficiency of the proton reconstruction and the trigger veto caused by the proton interaction with dead material is typically 98%, but goes down to 60% in the p x , p y region where the protons are expected to pass the RF shield present between the two RP stations.
TPC and TOF efficiencies were calculated as functions of particle p T , η and z vtx . The joint acceptance and reconstruction efficiency of a track in the TPC was measured separately for π + , π − , K + , K − , p andp using the single particle MC embedded into zero-bias trigger data taken simultaneously with physics triggers. TPC inefficiencies are caused by empty spaces between the sectors, fiducial cuts on the positions of space points, the influence of a high density of off-time tracks, the interactions of particles in dead material in front of the TPC, dead TPC modules, and natural decays of pions and kaons before or inside the TPC. The TPC efficiency increases with p T . The efficiency at the lowest transverse momentum used in the analysis is 70% for pions, 40% for kaons and 75% for protons. Above 1 GeV, the efficiency plateaus at 80% for pions, 70% for kaons and 85% for protons. The efficiency is roughly independent of particle charge, except that the efficiency for anti-protons is around 2% smaller than for protons. The TPC efficiencies depend on the track selection criteria. To check the sensitivity of the results on the track selection, the TPC efficiencies were calculated with looser and tighter matching criteria of the TPC tracks with a vertex and also with an increased (lowered) required number of associated hits 28 (20). These changes led to ±4% changes in the efficiencies.
The combined TOF acceptance, hit-reconstruction efficiency and matching efficiency with TPC tracks was measured separately for π + , π − , K + , K − , p andp using single particle MC embedded into zero-bias trigger data taken simultaneously with physics triggers. For high statistics exclusive π + π − production, the TOF efficiency was also measured using a data-driven tag-and-probe method where one of the pions matched with the TOF serves as a tag and the efficiency of the TOF was measured for the second pion. The differences between the data-driven and MC-based efficiencies were added to the MC-based efficiencies, assuming that they are independent of z vtx . The same procedure was applied to kaons and protons. The average p T -and η-dependent correction to the TOF efficiency is 3%. Finally, the TOF efficiency was measured directly by selecting in-time TPC tracks using an independent data sample with the HFT signal recorded. Reconstruction of TPC tracks containing hits in silicon layers of the HFT guarantees that tracks are in-time with the TOF hits. HFT matched tracks, however, have limited coverage of z vtx , with |z vtx | < 20 cm. The average additional correction is 1% for pions, 3% for kaons and 2% for protons.
A small fraction of signal events are rejected by the p miss T < 75 MeV requirement. The leakage is caused by the finite resolution of particle momenta which, in the case of a particle measured in the TPC, depends on the momentum. The efficiency of this cut was measured as a function of central particle momenta using a fast phase space MC simulation. This efficiency is 97% if both tracks have transverse momentum less than 0.4 GeV, and decreases to 89% for tracks with transverse momentum above 1.5 GeV.
The PID efficiency, defined as the probability that the particle pair h + h − passes the relevant PID selection criteria given by Eqs. (6.2), (6.3) and (6.4), was calculated as a function of the central particles' momenta using the fast phase space MC simulation. The PID efficiency for π + π − pairs is almost 100% in the whole fiducial region. For K + K − and pp it is also close to 100% if the lower of the two transverse momenta in the pair, p min T , is less than 0.6 GeV for K + K − or less than 1 GeV for pp. At larger values of p min T , efficiencies for K + K − and pp identification decrease significantly due to the χ 2 dE/dx (π + π − ) > 9 requirement used to limit misidentified π + π − pairs in the K + K − and pp samples.
The tracking efficiencies provide corrections for true particles inside the fiducial volume to be reconstructed in the TPC or RP detectors. Additional corrections were applied to account for true particles inside the fiducial volume that are reconstructed outside this volume, and true particles outside the fiducial region that are reconstructed inside this volume. Such migration is caused by finite detector resolutions and the intrinsic smearing of the forward proton kinematics due to the RHIC angular beam divergence. It is also possible that the presence of tracks not associated with the true particle causes an incomplete exclusive event to pass all the selection cuts. Such fake tracks may come from interactions of true particles with material in the detector or from additional pile-up processes. Correction factors for migrations through the boundary of the fiducial region, and for fake tracks, were estimated from GenEx and single particle samples. Joint correction factors for the migrations are generally very small, but can be up to 5% close to the edges of the fiducial region for central particles and up to 30% for forward protons. Corrections for fake particles only weakly depend on transverse momentum and are below 2% for central particles and up to 5% for forward protons.

Systematic uncertainties
Several sources of possible systematic uncertainties have been considered in this analysis. The largest contributions to the systematic uncertainty arise from the modelling of the RP system and the beam-line elements, detector alignment and the embedding technique. The overall uncertainty on measurement efficiencies related to the RP system is typically 6%, but up to 30% for |t 1 + t 2 | > 0.3 GeV 2 . This uncertainty is derived from the difference between efficiencies obtained from the MC simulation and from the data-driven tag-andprobe method using elastic scattering data.
The uncertainties related to the TPC efficiency are dominated by modelling of the disturbing activity in the detector caused by the high density of off-time tracks. An uncertainty of 1% was estimated by studying the consistency of the corrections obtained from the embedding technique for different rates of off-time tracks. An additional 1% uncertainty arises from the extrapolation of the embedding result, obtained from only a subset of data sample, to the full sample. Finally, a 0.5% uncertainty related to the amount of inactive material between the primary vertex and the STAR TPC was estimated based on the comparison of rates of secondary vertices between data and simulation. In addition, we observed up to ±1.5% changes in the cross sections by applying looser and tighter TPC track selection criteria and correcting using the TPC efficiency obtained for a given set of selection cuts. We treat these deviations as an additional source of systematic uncertainty. A typical TPC-related total systematic uncertainty on the cross section is 4% for π + π − and pp and 6% − 7% for K + K − .
The TOF-related uncertainties were estimated as the difference between results obtained with simulation and data-driven tag-and-probe methods using exclusive events and those obtained with the direct method using an independent sample of HFT-tagged tracks. A typical total TOF-related systematic uncertainty on the cross section is 3% for π + π − , 5% for pp and 10% for K + K − .
The only sources of systematic uncertainties which may vary significantly as a function of m(h + h − ) are non-exclusive and pile-up backgrounds. The quality of the description of these two sources of backgrounds is investigated using MC samples of Central Diffraction and Minimum Bias embedded into zero-bias collision data, using control samples enhanced in the background. The normalisation of both background sources for π + π − CEP are tuned to match the data in the signal-free regions of the control samples. The control sample for pile-up background normalisation was obtained by imposing all the standard selection cuts, except the requirement of consistency of the z-position of the vertex obtained from time difference of the signals from the forward protons in the RPs with the TPC vertex, ∆z vtx . The distribution of ∆z vtx is shown in Fig. 6 (left). The normalisation of the pile-up background shown by the red histogram is tuned to describe the data in the signal-free region of large difference between estimates of z-vertex. In the control sample for non-exclusive background normalisation, only the p miss T < 75 MeV cut was not used. The resulting p miss  for CEP π + π − event candidates between data and MC after offline selection excluding the total transverse momentum cut, marked with dashed red lines and arrows. In addition to the signal channel (opposite-sign particles), the control background channel (same-sign particles) is also shown in the plots. Data are represented by black (opposite-sign) or red (same-sign) points, while stacked MC predictions are drawn as filled (opposite-sign) or hatched (same-sign) histograms of different colours. Vertical error bars represent statistical uncertainties and horizontal bars represent bin sizes.
distribution is shown in Fig. 6 (right). The normalisation of the non-exclusive background shown by the green histogram is tuned to describe the data in the signal-free region of large values of p miss T . After tuning, the description of the data in both control regions is very good. The good agreement between data and simulation shown in Fig. 6 (right) is achieved only by removing the final states consisting of π + + π − + neutrals (mainly π + + π − + π 0 + π 0 ) from the PYTHIA8 CD prediction. Removal of π + + π − + neutrals not only makes the shape of the p miss T distribution compatible with data, but also correctly predicts the ratio of same-sign to opposite-sign pairs outside the signal region. Otherwise, this ratio is underestimated by 50%. The background in the signal sample is estimated by summing the predictions from both samples after applying all the selections cuts. Contributions of the remaining background in the data sample are shown in Fig. 7 together with estimations of the background using the (nominal) fully data-driven method, shown by magenta circles. For ∆ϕ (Fig. 7 (left)), the background contributions from the nominal and alternative methods agree well. Both methods show that the pile-up background is mainly located close to ∆ϕ = 180 • , as expected for pile-up from elastic scattering events. For m(ππ) (Fig. 7 (right)), the background contributions using the nominal and alternative methods disagree by up to 50% depending on the mass region. The enhancement of the data-driven estimate over the MC prediction around the f 2 (1270) and ρ 0 mass regions may be caused by imperfect modelling of the resonant states in the hadronisation model used in the PYTHIA generator. This discrepancy was not used as a systematic uncertainty on the background estimation. The third possible method for background estimation is another data-driven method using the same-sign control sample normalised to the opposite-sign sample in the signal-free region of p miss T > 0.15 GeV. The unscaled same-sign contribution is shown by red circles in Figs. 6 and 7, while the distributions predicted by PYTHIA 8 CD are shown by the blue hatched histogram. This method, by definition, accounts only for the combinatorial phase space background and cannot describe possible contributions from non-exclusive resonance production. This is seen in Fig. 7 (right), where the lack of ρ 0 and f 2 (1270) contributions is clearly seen. This method was also not used to evaluate systematic uncertainty. The systematic uncertainty related to background subtraction was estimated by replacing the polynomial function describing the shape of the background distribution with a histogram template obtained from the same-sign events normalised to opposite-sign events in the signal-free region of p miss T > 0.15 GeV. Agreement was found at the level of 10%, which is used as the systematic uncertainty. This contributes up to 1% uncertainty on the cross section.
The relative luminosity at STAR is determined from the coincidence rate in ZDC detectors in both beam directions. Absolute calibration is given by a special Van der Meer [57] scan. For the precise measurement of the total and elastic cross sections [49], three dedicated Van Table 1. Typical fractional systematic uncertainties of the integrated fiducial cross sections for CEP of π + π − , K + K − and pp pairs, decomposed into their major components.
uncertainty on the luminosity measurement. The luminosity uncertainty was estimated to be 4%. To account for possible fill-by-fill dependence in the luminosity measurement, an additional 4% uncertainty was assigned to the luminosity. It was determined by comparing variations of the effective cross sections for elastic scattering relative to the measurement done solely based on data collected during the fill with the Van der Meer scans. The overall luminosity uncertainty of 6% was estimated by the quadratic sum of uncertainties from these two sources.
Other systematic uncertainties considered include those due to the vertex reconstruction efficiency, selection cuts, and the trigger efficiency. None of these produce uncertainties beyond the 2% level.
The total systematic uncertainties are obtained by adding contributions from all sources in quadrature. They amount typically to between 10% and 20%, except at the extremes of the measurement range in |t 1 + t 2 |, and are highly correlated between bins. Table 1 shows the systematic uncertainties, decomposed into their major components, of cross sections for CEP π + π − , K + K − and pp pairs integrated over the fiducial region.

Cross sections in fiducial region
All results presented in this subsection are obtained in the fiducial region defined in Section 7. In Figs. 8, 9 and 10, the differential fiducial cross sections related to the variables characterizing the centrally-produced hadron pairs are presented. Figure 8 shows the differential cross section for CEP of π + π − pairs as a function of the pair invariant mass. There are several features of the distribution which need to be pointed out. The deep hole observed in the measured differential cross section dσ/dm(π + π − ) in the mass region m(π + π − ) < 0.6 GeV is mainly due to the fiducial cuts. At larger invariant masses, resonance structures are seen in the data consistent with the f 0 (980) and f 2 (1270) mesons expected to be produced in the Pomeron-Pomeron fusion process. At even higher invariant masses, another resonance is observed at ∼ 2.2 GeV. The DiMe model roughly describes both the normalisation and the shape of the continuum production under the resonances, up to masses of about 1.9 GeV. In contrast, the GenEx model fails to describe the shape of the continuum production. The MBR model prediction generally follows the shape of DiMe and GenEx predictions at masses below 1 GeV, but falls less rapidly with mass above 1 GeV.  Figure 8. Differential cross section for CEP of π + π − pairs as a function of the invariant mass of the pair in the fiducial region explained in the plot. Data are shown as solid points with error bars representing the statistical uncertainties. The typical systematic uncertainties are shown as dark/light gray boxes (with/without luminosity uncertainty included, respectively), for only a few data points as they are almost fully correlated between neighboring bins. Predictions from three MC models, GenEx, DiMe and MBR, are shown as histograms.
Notable are sharp drops of the predicted cross section at 0.7 GeV and 1.65 GeV. The former has been identified as a result of near-threshold-enhanced production of π + + π − + neutrals (mainly 2π 0 ), which starts in PYTHIA8 around 0.7 GeV. It has already been demonstrated in Section 10 that such events are extensively overpopulated in PYTHIA8. The latter drop of the cross section at 1.65 GeV, present also in the prediction for K + K − , results from the fiducial cut on central particle pseudorapidities |η| < 0.7 and peculiar correlation between the invariant mass and pseudorapidity of the final state particles in PYTHIA8. In Fig. 9, the differential fiducial cross sections for CEP of K + K − and pp pairs are shown. The measured differential cross section, dσ/dm(K + K − ), shows significant enhancement in the f 2 (1525) mass region and a possible smaller resonant signal in the mass region of f 2 (1270). Both structures are expected to be produced in the Pomeron-Pomeron fusion process. The ratio of the cross sections for π + π − to K + K − production in the f 2 (1270) mass region is roughly 18, consistent with the PDG ratio of the f 2 (1270) branching fractions for its decays into π + π − and K + K − [58], assuming similar contributions from non-resonant production under the f 2 (1270) peaks and similar STAR acceptance. The DiMe and GenEx predictions roughly describe the non-resonant contribution to the data in the resonance region. The data are also consistent with the ratio of the non-resonant exclusive production of π + π − to K + K − pairs expected by GenEx and DiMe. In the case of the differential cross section dσ/dm(pp), only predictions from the MBR model are available and they overestimate the data by a factor of 8.  Figure 10 shows the differential cross sections for CEP of different particle species pairs as a function of the pair rapidity. The shapes of the measured distributions are generally well described by all the model predictions.
In Fig. 11, the differential fiducial cross sections related to the forward-scattered protons are presented. Figure 11 (top) shows the differential cross sections for CEP of different particle species pairs as a function of ∆ϕ. Strong suppression of the differential cross sec-  Figure 10. Differential cross sections for CEP of charged particle pairs π + π − (left), K + K − (middle) and pp (right) as a function of the pair rapidity measured in the fiducial region explained in the Sec. 6. Data are shown as solid points with error bars representing the statistical uncertainties. The typical systematic uncertainties are shown as gray boxes for only a few data points as they are almost fully correlated between neighboring bins. Predictions from three MC models, GenEx, DiMe and MBR, are shown as histograms. In the lower panels, the ratios between the MC predictions (scaled to the data for better shape comparison) and the data are shown.   Figure 11. Differential cross sections for CEP of charged particle pairs π + π − (left column), K + K − (middle column) and pp (right column) as a function of the difference of azimuthal angles of the forward-scattered protons (top) and of the sum of the squares of the four-momenta losses in the proton vertices (bottom) measured in the fiducial region explained in the Sec. 6. Data are shown as points with error bars representing the statistical uncertainties. The typical systematic uncertainties are shown as gray boxes for only a few data points as they are almost fully correlated between neighboring bins. Predictions from three MC models, GenEx, DiMe and MBR, are shown as histograms. In the lower panels the ratios between the MC predictions (scaled to data) and the data are shown.
tions close to 90 • is due to the fiducial cuts applied to the forward scattered protons. The shape of DiMe model prediction agrees with data for π + π − and K + K − . The model implemented in GenEx does not describe the data. The MBR model implemented in PYTHIA8 describes the data fairly well in shape for K + K − and pp. Figure 11 (bottom) shows the differential cross sections for CEP of different particle species pairs as a function of |t 1 + t 2 |. The shapes of the measured cross sections are strongly affected by the fiducial cuts applied to the forward-scattered protons. The shapes of the differential cross sections for both π + π − and K + K − pair production are better described by the DiMe and MBR models than by the GenEx model. For pp pair production, the MBR model predicts a steeper slope.
The STAR detector acceptance naturally splits the fiducial region into two ranges of ∆ϕ, which are differently sensitive to absorption effects. Figure 12 shows the differential    Figure 12. Differential cross sections for CEP of charged particle pairs π + π − (top), K + K − (middle) and pp (bottom) as a function of the invariant mass of the pair in two ∆ϕ regions, ∆ϕ < 90 • (left column) and ∆ϕ > 90 • (right column), measured in the fiducial region explained on the plots. Data are shown as solid points with error bars representing the statistical uncertainties. The typical systematic uncertainties are shown as gray boxes for only a few data points as they are almost fully correlated between neighboring bins. Predictions from three MC models, GenEx, DiMe and MBR, are shown as histograms.   Figure 13. Differential cross sections dσ/dm(π + π − ) for CEP of π + π − pairs in two | p 1,T − p 2,T | regions, | p 1,T − p 2,T | < 0.12 GeV (left) and | p 1,T − p 2,T | > 0.12 GeV (right), in the fiducial region and with ∆ϕ < 90 • . There is no significant difference between the two | p 1,T − p 2,T | regions. Data are shown as solid points with error bars representing the statistical uncertainties. The typical systematic uncertainties are shown as gray boxes for only a few data points as they are almost fully correlated between neighboring bins. Predictions from three MC models, GenEx, DiMe and MBR, are shown as histograms.
cross sections for CEP of different particle species pairs as a function of the pair invariant mass in two ∆ϕ regions: ∆ϕ < 90 • (left column) and ∆ϕ > 90 • (right column). Sharp drops in the measured cross sections at m(π + π − ) < 0.6 GeV and at m(K + K − ) < 1.3 GeV for the ∆ϕ > 90 • range are due to the fiducial cuts applied to the forward-scattered protons.
In the case of the cross section for CEP of π + π − pairs in the ∆ϕ < 90 • range, the peak around the f 2 (1270) resonance in data is significantly suppressed while the peak at f 0 (980), as well as possible resonances in the mass ranges 1.3 − 1.5 GeV and 2.2 − 2.3 GeV, is enhanced compared to the ∆ϕ > 90 • range. Such correlations, between resonances observed in the mass spectrum and in azimuthal angle between outgoing protons, indicate factorisation breaking between the two proton vertices. In the range ∆ϕ < 90 • , the DiMe model describes well both the normalisation and the shape of the mass spectrum at m(π + π − ) < 0.5 GeV. In the cross section for CEP of K + K − pairs, the data do not show any significant ∆ϕ asymmetry except for a possible widening of the peak at f 2 (1520) in the region ∆ϕ < 90 • . This widening may indicate an enhancement of additional resonances around 1.7 GeV in this configuration. In the cross section for CEP of pp pairs, the data do not show a significant ∆ϕ asymmetry except for a possible enhancement in the 2.2 − 2.4 GeV mass range for the ∆ϕ > 90 • region. Experimental observation of vertex factorisation breaking in the pp collisions motivated Close and Kirk in Ref. [59] to propose a method for filtering glueballs from their qq counterparts. The gg configurations were proposed to be enhanced in the limit | p 1,T − p 2,T | → 0. Such a configuration is already enhanced in the ∆ϕ < 90 • region. To further enhance a possible gg configuration, the data are studied as a function of | p 1,T − p 2,T | in the ∆ϕ < 90 • region. Figure 13 shows the differential cross section for CEP of π + π − as a function of the pair invariant mass separately in two | p 1,T − p 2,T | regions: | p 1,T − p 2,T | < 0.12 GeV(left)  Figure 14. Differential cross sections for CEP of charged particle pairs π + π − (left column), K + K − (middle column) and pp (right column) as a function of cos θ CS (top) and of φ CS (bottom), measured in the fiducial region explained in the Sec. 6. Data are shown as solid points with error bars representing the statistical uncertainties. The typical systematic uncertainties are shown as gray boxes for only a few data points as they are almost fully correlated between neighboring bins. Predictions from three MC models, GenEx, DiMe and MBR, are shown as histograms. In the lower panels the ratios of the MC predictions scaled to data and the data are shown. and | p 1,T − p 2,T | > 0.12 GeV(right). The data do not show any changes in the shape of the π + π − mass spectrum for the two ranges of | p 1,T − p 2,T | after filtering events with ∆ϕ < 90 • .
We have also studied the angular distributions of the charged particles produced in the final state, which may help to constrain the underlying reaction mechanism. This can be done in various reference frames. However, for an easy comparison with theoretical predictions, we use here the Collins-Soper [56] reference frame which is also used, e.g., in Ref. [60]. Figure 14 (top) shows the differential cross sections for CEP of different particle species pairs as a function of cos θ CS . In general, the model predictions are narrower than the data for all particle species pairs. The only exception is the DiMe prediction for π + π − production, which fits the data much better than other models. Figure 14 (bottom) shows the differential cross sections for CEP of different particle species pairs as a function of φ CS . None of the models is able to describe the data. The double peak structure observed in the data is due to the STAR TPC acceptance.
High statistics of the two-pion sample allow to study the CEP of π + π − pairs in greater detail. Figure 15 shows the differential cross sections for CEP of π + π − pairs as a function  Figure 15. Differential cross sections for CEP of π + π − pairs as a function of the rapidity of the pair (left column), the difference in azimuthal angles of the forward-scattered protons (middle column) and the sum of the squares of the four-momentum losses in the proton vertices (right column) measured in the fiducial region explained in the Sec. 6, separately for three ranges of the π + π − pair invariant mass: m < 1 GeV (top), 1 GeV < m < 1.5 GeV (middle) and m > 1.5 GeV (bottom). Data are shown as points with error bars representing the statistical uncertainties. The typical systematic uncertainties are shown as gray boxes for only a few data points as they are almost fully correlated between neighboring bins. Predictions from three MC models, GenEx, DiMe and MBR, are shown as histograms. In the lower panels, the ratios between the MC predictions (scaled to data) and the data are shown. of the pair rapidity (left column), ∆ϕ (middle column) and |t 1 + t 2 | (right column) in three characteristic ranges of the invariant mass of the pair: m(π + π − ) < 1 GeV (mainly nonresonant production), 1 GeV < m(π + π − ) < 1.5 GeV (f 2 (1270) mass range) and m(π + π − ) > 1.5 GeV (higher invariant masses). In the case of the cross section, dσ/dy, all the models agree with the shape of the data in all three mass ranges except for the GenEx and DiMe predictions in the highest mass range, where the predictions are narrower than the data. Strong suppression of the fiducial cross section close to ∆ϕ = 90 • is due to the STAR RP acceptance, while the asymmetry of ∆ϕ = 0 • vs. ∆ϕ = 180 • in the lowest mass region is due to the STAR TPC acceptance. The DiMe model agrees with data only in the lowest mass range. The model implemented in GenEx does not describe the data in any of the three mass regions. Both DiMe and GenEx show similar shapes in the ∆ϕ distribution except in the lowest mass region. The MBR model predicts symmetric ∆ϕ distributions in all mass ranges, which is not supported by the data except in the highest mass region. The slope of the cross section as a function of |t 1 + t 2 | is less steep in the f 2 (1270) mass range compared to other mass ranges. A comparison between the model predictions of the |t 1 + t 2 | distributions and the data does not show a significant mass dependence. The best description is given by the MBR model, and the worst by GenEx. Figure 16 shows the differential cross sections for CEP of π + π − pairs as a function of cos θ CS (top row) and φ CS (bottom row) in three characteristic ranges of the invariant mass of the pair: m(π + π − ) < 1 GeV, 1 GeV < m(π + π − ) < 1.5 GeV and m(π + π − ) > 1.5 GeV. To help in interpreting the data and in understanding the STAR acceptance effects, the data are compared with expectations (like angular distributions of pions in the π + π − rest frame) from models with pure S 0 and D 0 waves. The S 0 wave predicts a uniform distribution of polar angle φ, in contrast to the D 0 wave. The angular distributions are generated in the most natural Gottfried-Jackson frame [61] with the Pomeron-Pomeron direction taken as the z-axis. The transformation to the Collins-Soper frame changes the angular distributions for the D 0 wave but not for the S 0 wave. Therefore, the shape of the φ CS distribution for S 0 wave, after applying fiducial cuts, represents also the φ CS shape of the STAR acceptance. The S 0 and D 0 predictions are normalised to data. The double-peak structure observed in the φ CS distribution in the lowest mass region, where data are reasonably well described by S 0 prediction, is due to the STAR acceptance. In contrast, at higher masses, where prediction from the S 0 wave model is flat, the double-peak structure does not come from the STAR acceptance. Both cos θ CS and φ CS in the lowest mass region agree very well with the S 0 wave suggesting that this mass region is dominated by spin-0 contribution. At higher masses, pure S 0 or D 0 waves are not able to describe the data.
In the case of the differential cross section dσ/d cos θ CS , the DiMe predictions fit the data only in the lowest mass region. In contrast, the MBR predictions fail to describe the shape of the cos θ CS distribution in this mass range only. The GenEx prediction does not describe the data in any mass range. In the case of the differential cross section dσ/dφ CS , in the lowest mass region only GenEx predicts the shape of the φ CS distribution. The DiMe prediction fits the data well in the middle mass range. Both GenEx and DiMe predictions describe the shape of the φ CS distribution fairly well in the highest mass region.
The cross sections, integrated over the full fiducial range of the analysis, are shown in   Figure 16. Differential cross sections for CEP of π + π − pairs as a function of cos θ CS (top) and of φ CS (bottom) measured in the fiducial region is shown in three ranges of the π + π − pair invariant mass: m < 1 GeV (left column), 1 GeV < m < 1.5 GeV (middle column) and m > 1.5 GeV (right column). Data are shown as solid points with error bars representing their statistical uncertainties. The typical systematic uncertainties are shown as gray boxes for only a few data points as they are almost fully correlated between neighboring bins. Predictions from three MC models, GenEx, DiMe and MBR, as well as from pure S 0 and D 0 waves are shown as histograms. In the lower panels the ratios between the MC predictions (scaled to data) and the data are shown.  Table 2. Integrated fiducial cross sections with statistical and systematic uncertainties for CEP of π + π − , K + K − and pp pairs in two ranges of azimuthal angle difference, ∆ϕ, between the two forward-scattered protons. Table 2 in two ranges of ∆ϕ. The largest contribution to the uncertainty of π + π − production arises from the luminosity measurement. For K + K − production, the largest con-tribution to the uncertainty arises from the TOF efficiency. In case of pp production, the uncertainty is dominated by statistical fluctuations.
11.2 Extrapolated π + π − differential cross sections dσ/dm and d 2 σ/dt 1 dt 2 Invariant mass distributions in the fiducial region of the measurement cannot be directly used to extract yields of possible resonances without extrapolation to the full kinematic region of the central pion pair, given by p T → 0 and |η| → ∞ (full solid angle in the central system rest frame). Extrapolation to an unmeasured region is always model dependent. In this section, we present the cross section corrected to the full phase space using a flat angular approximation which distributes scalar decays uniformly over the solid angle in the rest frame of the central system. This choice is supported by the generally good description of the pion angular distribution by the S 0 wave distribution shown in Fig. 16, and by the expected dominant production of 0-spin states. However, other scenarios are also considered. To limit the corrections, the measurement is restricted to |y(π + π − )| < 0.4. This keeps scalar decays uniform, by Lorentz invariance. In the correction calculation, the factorisation of the phase space of the central system and forward protons is assumed. For the forward protons' phase space, a uniform distribution of azimuthal angles is assumed, while polar angles are generated according to an exponential t distribution with t-slope of 6 GeV −2 . The measurement is extrapolated from the part of the fiducial region given by Eq. (6.1), covering 0.05 GeV 2 ≤ −t 1 , −t 2 ≤ 0.16 GeV 2 , to the Lorentz-invariant phase space region defined by the same t interval and full azimuthal angle of forward-protons. The measurement is further restricted to two ranges of ∆ϕ, ∆ϕ < 45 • and ∆ϕ > 135 • , which reduces the extrapolation correction and the systematic uncertainties related to their modelling. A minimal model of the π + π − invariant mass spectrum was fitted to the extrapolated differential cross section. In this model, we assume contributions from direct pair production and only three resonances in the mass range of 0.6−1.7 GeV: f 0 (980), f 2 (1270) and f 0 (1500). The total amplitude for the exclusive π + π − production is then given by Thus all states are added coherently and can interfere with each other. The amplitude for continuum production is chosen to be real, while the amplitudes for the production cross sections for resonances in the π + π − channel are allowed to have non-zero phase shifts, φ. The shape of the continuum amplitude in the fitted mass range is assumed to have the form with the break-up momentum q(m) = 1 2 m 2 − 4m 2 π . This continuum effectively includes the production of other wide resonant states, e.g. f 0 (500), with a phase of amplitude slowly varying within the fit range, as described below in the discussion of the fit result. For the f 2 (1270) and f 0 (1500) resonances, we use the relativistic Breit-Wigner form of the production amplitude with mass-dependent width: A factor I is introduced to provide proper normalisation, +∞ 2mπ dm|R BW | 2 = 1. As a result, the total cross sections, σ f 0 (980) , σ f 2 (1270) and σ f 0 (1500) , for resonance production are directly obtained from the fit. The centrifugal effects are accounted for in Eq. (11.3) through the Blatt-Weisskopf barrier factors, B J [62,63], with the empirical interaction radius, R, set to 1 fm and q 0 = q(M ). J = 2 and J = 0 are used for f 2 (1270) and f 0 (1500), respectively.
The f 0 (980) meson requires a different treatment due to the large branching ratio to the KK channel, which opens in the vicinity of the mass peak. This changes the resonance shape and is accounted for in the parameterisation of the amplitude via the Flatté formula [64]: . (11.4) The partial widths, Γ j (j = π, K), are described by the product of the coupling parameter g j and the break-up momentum q j for particle j: is the partial width in the π + π − channel at the resonance mass. In the fit, the ratio g K /g π is fixed to 4.21, the value well-constrained experimentally through the measurement of J/ψ decays into φ mesons and π + π − /K + K − pairs [65].
To fit model parameters to the data, the binning in the invariant mass distribution must be adjusted to the expected structures in the cross section. This requires narrower binning than for the measurement of the fiducial cross section and taking account of the impact of detector resolution. The squared amplitude from Eq. (11.1), |A| 2 , is convoluted for the purpose of the fit with the normal distribution, N (0, σ res ), representing the finite measurement resolution of the invariant mass of the pion pair. The resolution parameter, σ res (m), is provided to the fitting algorithm; it is set to grow linearly with increasing invariant mass according to MC simulation of the STAR TPC detector. The m(π + π − ) resolution at the lower and upper limit of the fit range is equal to 4 MeV and 13 MeV, respectively. The final form of the function fitted to the extrapolated dσ/dm(π + π − ) distribution is given by the convolution of the total amplitude squared with the normal distribution: The fitting is performed using the Minuit2 toolkit [66] within the ROOT analysis software [67]. The standardly-defined χ 2 is minimised simultaneously in two ∆ϕ ranges. For each of the two f 0 resonances the fitted values of mass and width in the two ∆ϕ subsets are forced to be equal, while the phases and absolute values of amplitudes of all resonances are left independent. The mass and width of the f 2 (1270) resonance is fixed to the well-known Particle Data Group values [58].
The experimental systematic uncertainties of the model parameters are estimated through the independent fits to the extrapolated dσ/dm(π + π − ) distributions with each of the systematic variations described in Sec. 10 applied. In addition to this, we take into consideration the sensitivity of the fit result to the modelling of the extrapolation to the full kinematic region. We check the effect of the extrapolation to the full solid angle in the π + π − rest frame, assuming a smooth transition from the angular distributions for pure S 0wave (up to 1 GeV) to the angular distributions for pure D 0 -wave (starting at 1.2 GeV). We also check the effect of using the extrapolation calculated with predictions from the DiMe and GenEx generators, both for the central state and for the forward-scattered protons. We also check the result of the fit with the ratio g K /g π varied within its uncertainties. The systematic uncertainty on a parameter is separated into two parts. The first one is related to the experimental uncertainties and is calculated as the quadratic sum of the differences between the nominal fit result and the result of the fit to dσ/dm(π + π − ) with each systematic effect. The second part is related to the extrapolation and is quoted as the largest deviation from the nominal result.
The extrapolated cross sections are shown in Fig. 17, together with the result of the fit described above. The model parameters providing the minimum χ 2 are listed in Tab. 3. The fit, with a total of 20 free parameters, gives χ 2 /ndf = 175/180 which shows that the data and the model are in excellent agreement in the fit region. Alternative extrapolation models show a similar fit quality, although some parameters change significantly as can be noted from the model-related uncertainties in Tab. 3. The fitted model shows a small deviation from the extrapolated data around 1.37 GeV. This might result from the presence of the f 0 (1370), however the inclusion of the f 0 (1370) is not necessary to describe the data. The cross section for f 0 (1500) production differs from zero by 5 and 2 standard deviations in the ∆ϕ < 45 • and ∆ϕ > 135 • regions, respectively. Removing the f 0 (1500) from the fit makes the χ 2 /ndf change to 352/186, a 7.0 standard-deviation effect. From this, we infer that the shape of dσ/dm(π + π − ) around 1.4 − 1.6 GeV, the high-mass part of the f 2 (1270) region, is primarily determined by the f 0 (1500) interfering with π + π − continuum.
Since the masses and widths of the f 0 (980) and the f 0 (1500) are free parameters, we can compare their fitted values with the PDG data [58]. In the case of the f 0 (980), the mass and width are found to be M f 0 (980) = 956 ± 7(stat.) ± 1(syst.) +4 −6 (mod.) MeV and Γ 0,f 0 (980) = 163 ± 26(stat.) ± 3(syst.) +17 −20 (mod.) MeV, respectively. These values differ from the PDG estimates of mass (990±20 MeV) and width (from 10 MeV to 100 MeV). However, the PDG emphasises a strong dependence of the resonance parameters on the model of the amplitude. Some measurements listed in Ref. [58] are in reasonable agreement with our measured numbers. In addition to this, the mass and width of the f 0 (980) resulting from the fit with the Breit-Wigner form (instead of the Flatté form) of the amplitude gives a result M f 0 (980) = 974 ± 1(stat.) ± 1(syst.) MeV and Γ 0,f 0 (980) = 65 ± 3(stat.) ± 1(syst.) MeV (albeit with a notably worse χ 2 /ndf of 226/180 providing evidence for a significant branching   Figure 17. Differential cross section dσ/dm(π + π − ) extrapolated from the fiducial region to the Lorentz-invariant phase space given by the central-state rapidity, |y(π + π − )| < 0.4, and squared four-momentum transferred in forward proton vertices, 0.05 GeV 2 < −t 1 , −t 2 < 0.16 GeV 2 . The left and right columns show the cross sections for ∆ϕ < 45 • and ∆ϕ > 135 • , respectively. The data are shown as black points with error bars representing statistical uncertainties. The result of the fit, F(m), is drawn with a solid red line. The squared amplitudes for the continuum and resonance production are drawn with lines of different colours, as explained in the legend. The most significant interference terms are plotted in the middle panels, while the relative differences between each data point and the fitted model is shown in the bottom panels.
fraction for the decay into KK, which needs to be accounted for in the resonance parameterisation). These values are in excellent agreement with PDG estimates and f 0 (980) parameters from other measurements that assume a Breit-Wigner resonance shape [58]. For the f 0 (1500), we obtain from the fit M f 0 (1500) = 1469 ± 9(stat.) ± 1(syst.) ± 2(mod.) MeV and Γ 0,f 0 (1500) = 89 ± 14(stat.) ± 2(syst.) +4 −3 (mod.) MeV. These values also deviate from the PDG averages 1505 ± 6 MeV for the mass and 109 ± 7 MeV for the width. However, numerous measurements on f 0 (1500) referenced in the PDG (and not used for the average calculation) report masses below 1500 MeV and widths below 100 MeV that are consistent with our result.
We have tested the possibility of the existence of an additional resonance produced in the mass range 1.2 − 1.5 GeV. With an f 0 -like component added to the model in Eq. (11.1), the best fit is achieved for M f 0 = 1372 ± 13(stat.) MeV and Γ 0,f 0 = 44 ± 24(stat.) MeV. In that case, the χ 2 /ndf is equal to 158/174 (p-value: 0.8), compared to the nominal value 175/180 (p-value: 0.6). The cross section dσ/dm(π + π − ) around 1. 37 Table 3. Results of the fit described in the text in two ranges of azimuthal angle difference ∆ϕ between forward-scattered protons. Statistical, systematic and model uncertainties are provided for each parameter, in that order.
is better described than with the nominal fit. Other parameters in the fit change slightly but remain compatible with their original values. The fitted content of the additional f 0 resonance is several times lower than the extracted yield of f 0 (1500) for ∆ϕ < 45 • , while for ∆ϕ > 135 • it is consistent with 0. The mass agrees with that of the f 0 (1370) resonance, however the measured width is much narrower than PDG estimates of about 200−500 MeV. We calculated the ratios of the total cross sections σ f 0 (980) /σ f 2 (1270) and σ f 0 (1500) /σ f 2 (1270) in two ∆ϕ regions, as well as the ratio σ(∆ϕ < 45 • )/σ(∆ϕ > 135 • ) for all resonances, as listed in Tab. 4. In the ratios, many of the systematic uncertainties cancelled out. We observe a significant dependence of the resonance production cross sections on the azimuthal separation of the forward-scattered protons. The two scalar mesons, f 0 (980) and f 0 (1500), are produced predominantly at ∆ϕ < 45 • , whereas the tensor meson f 2 (1270) is produced predominantly at ∆ϕ > 135 • . This ∆ϕ dependence is consistent with the observation   Figure 18. Comparison of the continuum production cross section obtained from the fit to the extrapolated dσ/dm(π + π − ) (dotted black line) with predictions from the continuum models. Scaled GenEx predictions are shown with the blue histogram. Three predictions from DiMe representing different meson form factors are shown as coloured bands, each spanning between minimum and maximum predicted cross sections from four available absorption models. Left and right panels show cross sections for ∆ϕ < 45 • and ∆ϕ > 135 • , respectively. made by the WA102 Collaboration [68].
The differential cross section for the π + π − continuum production, dσ cont /dm(π + π − ) = A 2 cont f 2 cont (m), is extracted from the fit to the extrapolated dσ/dm(π + π − ) distribution. It is then compared with expectations from GenEx and DiMe models using all three forms of meson form factors (Λ 2 off = 1.0 GeV 2 , b = 2 GeV −1 and a 2 = 0.5 GeV 2 , and a 0 = 1.0 GeV 2 ) and four models of absorption [19]. All models predict a shape for dσ cont /dm(π + π − ) that is different from the assumed form, f 2 cont , as shown in Fig. 18. The shape of the continuum predicted by the models can be achieved by changing the factor q/m in Eq. (11.2) to ( q/m) P , with the parameter P taking values between 2 and 10. Using such a modified continuum amplitude in the fit one obtains a P parameter consistent with 1 in both ∆ϕ ranges, and the remaining parameters are consistent with the results of the nominal fit. The deviation of the fitted dσ cont /dm(π + π − ) from all the model predictions is most evident at the lower edge of the fit range. A possible explanation of this observation is the presence of the f 0 (500) state, expected in DPE, and the photo-produced ρ 0 vector meson, generally suppressed within the kinematic region of this measurement. One should therefore treat the continuum obtained from the fit as an effective description of the coherent sum of the continuum and other states not explicitly included in the Eq. (11.1).
For these values, we do not separate modelling uncertainties since they are generally much smaller than experimental uncertainties. This is a consequence of the uniform ϕ distribution in all the models and the rather weak dependence of the cross section on ∆ϕ within the measured ranges. Such approximations are well-founded and in good agreement with the data. Variations of the slope, β, with m(π + π − ) and ∆ϕ can give important constraints for model developers. For example, it was pointed out in Ref. [32] that the f 2 (1270) cross section may show an enhancement when t 1 → 0 and t 2 → 0 for some couplings, while for others a suppression is expected. This enhancement or suppression results in a larger or smaller value of β, respectively.

Summary
We have studied the CEP of charged particle pairs (π + π − , K + K − and pp) in events with forward protons tagged in the RP detectors, using a sample of 14.2 pb −1 data collected with the STAR detector at RHIC in proton-proton collisions at √ s = 200 GeV. The centreof-mass energy of pp collisions in the present measurement is three times larger than the previous highest-energy measurement of the DPE process with forward-proton tagging performed at CERN at the ISR in the AFS experiment [25]. The uncertainty of the absolute normalisation of the STAR measurement is a factor of four better compared to measurements at the ISR, giving much stronger constraints for phenomenological models. The fits to the extrapolated cross sections as a function of the π + π − invariant mass of the minimal model, including the f 0 (980), f 2 (1270) and f 0 (1500) resonances and the non-resonant contribution, provide for the first time measurements of the relative phases between all the production amplitudes. In the glueball sector, there is no evidence except for a small enhancement around 1.7 GeV for production of f 0 (1710) decaying to a pair of kaons, and only weak evidence for f 0 (1370) production decaying to a pair of pions. However, we measured the production of a scalar meson which could have a gluonic contribution mixed in (f 0 (1500)). Measured masses and widths of the f 0 (980) and f 0 (1500), together with extensive studies of the non-resonant "background", may provide constraints for better understanding the role of the gluonic component in the family of scalar mesons.

Cross sections in fiducial region
The mass spectrum of the π + π − pairs shows an order-of-magnitude drop at 1.0 GeV, a clear peak around 1.3 GeV and possible further structures at higher masses. This is consistent with expectations for the DPE process and with experiments at lower energies investigating the production of f 0 (980), f 2 (1270) and also higher-mass resonances. The mass spectrum of K + K − shows a clear peak above 1.5 GeV and possible enhancement around 1.3 GeV. This is also consistent with expectations for the DPE process and with experiments at lower energies investigating the production of the f 2 (1270) and f 2 (1520), and possibly also f 0 (1710) resonances. The ratio of the measured cross sections for π + π − to K + K − production in the f 2 (1270) mass region is roughly 18, assuming a similar contribution from non-resonant production under the f 2 (1270) peak and similar STAR acceptance. This is consistent with the PDG ratio of f 2 (1270) branching fractions for decays into π + π − and K + K − pairs.
The mass spectrum of π + π − shows that, for ∆ϕ < 90 • , the peak around the f 2 (1270) resonance is significantly suppressed while the peak at f 0 (980), as well as possible resonances in the mass range 1.3 − 1.5 GeV, are enhanced compared to the region of ∆ϕ > 90 • . Such a correlation between resonances seen in the mass spectrum and the azimuthal angle between outgoing forward protons indicates factorisation breaking between the two proton vertices, reported for the first time by the WA91 experiment [26]. The present data do not show significant changes in the shape of the π + π − mass spectrum as a function of | p 1,T − p 2,T | after filtering events with ∆ϕ < 90 • . The mass spectra of the K + K − and pp pairs also show some indications of factorisation breaking.
The |t 1 + t 2 | distribution in the case of the π + π − pair production in the f 2 (1270) mass region is less steep compared to other mass ranges, suggesting that the t-slope of f 2 (1270) production is smaller than that of other states.
The angular distributions of pions in the π + π − rest frame indicate that, for invariant masses below 1 GeV, the data are dominated by spin-0 contributions. In the higher-mass region, the data show significant contributions from higher-spin states.
The measured cross sections are compared to phenomenological predictions from DPE models implemented in the form of the Monte Carlo generators GenEx, DiMe and PYTHIA8 (MBR model). The cross sections for CEP of π + π − and K + K − pairs are significantly above the GenEx and DiMe predictions. This is expected, as neither of these models includes contributions from resonant production. The shapes of distributions other than the mass spectra show generally good agreement between the data and predictions, especially for DiMe. GenEx predicts shapes slightly different from the data for the ∆ϕ, cos(θ CS ) and y(π + π − ) distributions. This might be attributed to absorption effects, which are treated fully differentially in DiMe and only on average in GenEx.
The shapes of several distributions for pp production are reasonably well described by PYTHIA8 (MBR model) but the prediction overestimates the data by a factor 8. The limited statistics do not allow any significant conclusions about the expected resonances below m(pp) ≈ 2.5 GeV.
The phenomenological interpretation of the data requires improvements of the DPE models to consistently include the continuum and resonance-production mechanisms, and the interference between the two, as well as absorption and rescattering effects.

Extrapolated differential cross sections
The fiducial measurement of π + π − production was extrapolated to the Lorentz-invariant phase space given by |y(π + π − )| < 0.4 and 0.05 GeV 2 < −t 1 , −t 2 < 0.16 GeV 2 in two ranges of ∆ϕ < 45 • and ∆ϕ > 135 • . This allows us to fit the extrapolated differential cross section with a minimal model of the π + π − invariant mass spectrum consisting of f 0 (980), f 2 (1270) and f 0 (1500) resonances and direct non-resonant π + π − production. The masses and widths of the f 0 (980) and f 0 (1500) resonances obtained from the fit are in good agreement with the PDG values. The two scalar mesons, f 0 (980) and f 0 (1500), are predominantly produced at ∆ϕ < 45 • , whereas the tensor meson f 2 (1270) is predominantly produced at ∆ϕ > 135 • . We observed weak evidence for an additional resonant state with a mass of 1372 ± 13(stat.) MeV and a width of 44 ± 24(stat.) MeV. This mass agrees with that of the f 0 (1370), but the width is much narrower than the PDG value of 200 − 500 MeV. The extrapolated cross sections of continuum production within the mass range 0.6 GeV < m < 1.7 GeV show an expected ∆ϕ asymmetry caused by absorption. Fits to an exponential function of the form ∝ exp [βt 1 ] exp [βt 2 ] were performed on the differential cross sections d 2 σ/dt 1 dt 2 to extract the slope of the −t distributions. Variations in the slope with m(π + π − ) and ∆ϕ can give important constraints for construction of phenomenological models of CEP of π + π − pairs.