Charged particle pseudorapidity distributions measured with the STAR EPD

In 2018, in preparation for the Beam Energy Scan II, the STAR detector was upgraded with the Event Plane Detector (EPD). The instrument enhanced STAR's capabilities in centrality determination for fluctuation measurements, event plane resolution for flow measurements and in triggering overall. Due to its fine radial granularity, it can also be utilized to measure pseudorapidity distributions of the produced charged primary particles, in EPD's pseudorapidity coverage of $2.15<|\eta|<5.09$. As such a measurement cannot be done directly, the response of the detector to the primary particles has to be understood well. The detector response matrix was determined via Monte Carlo simulations and corrected charged particle pseudorapidity distributions were obtained in Au+Au collisions at center of mass collision energies $\sqrt{s_{_{NN}}}$ = 19.6 and 27.0 GeV using an iterative unfolding procedure. Several systematic checks of the method were also done.

measurements, independent centrality determination for fluctuation measurements, and using it as a trigger in high luminosity environment during the BES-II program.
The EPD is a completely new subdetector that was supposed to improve the event plane resolution: for example, by about a factor of 2 in Au+Au collisions at √ s NN = 19.6 GeV [3]. Its predecessor (in event plane determination), the Beam-Beam Counter (BBC) has much less fine granularity than the EPD: only 36 tiles, with the 18 inner smaller tiles used -compared to the 372 tiles of the EPD [2]. It also has smaller acceptance of 3.3 < |η| < 5.0 in pseudorapidity [4]. The detector consists of two "wheels" on either (West and East) side of the STAR detector system, installed ±375 cm from the nominal interaction point (the detector's center). Each wheel consists of 12 "supersectors" covering φ = 30 • in azimuthal angle, each further segmented to 31 "tiles", thus giving 16 radial segments so-called "rings" 2 covering relatively large forward pseudorapidity range of 2.15 < |η| < 5.09 (or, range of 0.7 • < θ < 13.5 • angle to particle beam axis). Each supersector is connected to a bundle of 31 optical cables that transport light to high-efficiency silicone photo-multipliers (SiPM). The signals are then sent to the digital data acquisition systems [2].
Each tile registers hits, mostly Minimum Ionizing Particles (MIPs). Assuming that the probability distribution of the measured signal of a single hit can be described by a Landau distribution, the presence of multiple hits will result in a convolution of multiple Landau distributions. The measured Analog Digital Count (ADC) distributions were fitted with a multi-MIP Landau function, shown in Fig. 1. The different Landau distributions corresponding to the ADC contribution caused by n number of MIPs were convolved with different convolution weights (n-MIP weight).
The conclusion drawn was that convolving with less than 5 n-MIP weights are adequate to achieve a good fit, as the contribution of the 5-MIP weight was already zero within uncertainties -under the asssumption that the MIP weights were Poisson-distributed which was validated during data analysis. In view of this result, the systematic uncertainty contribution from this source -that is, fitting only up to 5 n-MIP weights -can be considered negligible.

Charged particle pseudorapidity measurement with the EPD
The aim is to measure the angle θ between the three-momentum p of the particle and the beam. Instead, a more convenient 3 quantity, the pseudorapidity η is used, which is defined as: where p z is the z component of the momentum, and the z direction is chosen to coincide with the direction of the beam [5]. Beyond the event plane determination, the EPD's fine radial granularity allows for pseudorapidity measurements to be performed. The raw EPD hit numbers could be used to calculate the pseudorapidity distribution of charged particles (dN ch /dη) by using the corresponding η value of the given ring. 2 The rings are numbered from 1 to 32 in the following manner: the innermost East EPD ring is the #1 which follows outerwards until #16; then, the #17 continues on the West EPD side's outermost ring until #32 being the innermost one. 3 In the ultrarelativistic limit, it approaches to rapidity (in c = 1 unit system, c being the speed of light): η ≈ y ≡ 1 2 ln E+pz E−pz , with E being the energy of the particle. However, this also includes the secondary particles that do not originate from the primary vertex. As the EPD is preceded by the rest of the detector system and is relatively far from the interaction point, multiple factors distort ("blur") the measured distribution.
The factors assumed to cause the most significant distortion effect are as follows. First of all, charged primary particles scatter in detector material (or in rare cases with each other), creating secondary particles contributing to dN ch /dη significantly. This is demonstrated in Fig. 2a, where the vertices (origins) of particles hitting the EPD in a detector material simulation are depicted. Secondly, neutral primary particles contribute through decays (e.g. a neutral Λ baryon decaying into proton and pion). In Fig. 2b it is clearly demonstrated that this contribution is non-negligible (based on the same simulation as mentioned above).

From raw EPD data to pseudorapidity distribution [dN/dη]
Using the previously mentioned multi-MIP Landau fit, one can extract the number of EPD hits for each ring; denoted as N (i Ring ) in the i-th ring. Given the underlying pseudorapidity distribution of the primary particles (dN/dη), assuming linear dependence from the dN/dη, the number of hits in a given ring can be calculated formally as a convolution: where R denotes the response matrix, which encodes response of the detector, i.e. connects a detector-level distribution with the true distribution to be measured. In this analysis, it contains the number of hits in the given ring number distribution's bin, originating from a particle at given η pseudorapidity distribution's bin.
No probabilistic consideration guarantees this matrix to be invertible, therefore a simple (or even a regularized) matrix inversion might not be an option even if the exact form of R would be known. Instead, a method called bayesian iterative unfolding [8] ("deblurring") is used.
Using this approach, the R needs to be extracted from simulations that are as close to the real system as possible. Using a complex event generator, a list of primary particles is obtained, along with a list of EPD hits -preferably all linked to primary tracks causing them. In this analysis, the events were generated using the STAR's HIJING Monte Carlo event generator combined with Geant4 to simulate the precise geometry of the EPD. In the following, the abbreviation MC will indicate data from these simulations. Such a response matrix can be seen in Fig. 3.
It should be noted that no (light) ion fragments can be simulated in HIJING, which are, in reality, inevitable with heavy-ion collisions. However, this shortfall should not change the results significantly, according to PHOBOS results [9]: the contribution from light ion fragments causes at least an order of magnitude smaller contribution to dN/dη than the results in this analysis (see Sec. 4).
In the following step, the unfolding technique is utilized to determine an uncorrected dN /dη. The software used for this purpose is the RooUnfold [10] framework, implemented in C++, running within the ROOT environment [11]. The package itself defines classes for the different unfolding algorithms -amongst others, the bayesian iterative unfolding.
The response matrix class of the software includes functions for populating the response matrix 4 as well as for managing the background (missed hits from real primaries and hits resulting from other sources 5 ).
During the unfolding, one can choose to propagate the statistical uncertainty in different ways; in this case, the most appropriate method should be propagating the (mostly badly conditioned, thus non-invertible) covariance matrix [8].
The resulting EPD ring distribution 6 needs to be corrected for the multiple counting (efficiency, ), explained as follows. The unfolding procedure results in one unfolded track for each individual EPD hit. However, it should be noted that one primary track can cause multiple hits. This effect needs to be corrected for -either via a bin-by-bin correction calculated 4 Fill(x measured , x truth ); naturally, "measured" and "truth" here stand for the training datasets obtained from MC (simulation). 5 Miss(x truth ) and Fake(x measured ) 6 Caused by both primary and secondary particles. Figure 3: Heatmap visualization of the R response matrix, connecting bins containing numbers of EPD ring hits (caused by either primary or secondary particles) with bins corresponding to primary particles at given η pseudorapidity. The left side corresponds to East EPD wheel, the right side to West EPD wheel. It is worth noting that many primaries create hits even in the opposite side EPD via secondaries, as seen in upper left and bottom right quarters.
from MC data (via a Number of hits from 1 primary(η) distribution), or by weighing the values filled in response matrix such that it could compensate for the multiple counts during the unfolding. In this analysis, the first method was used.

Extracting charged particle pseudorapidity distribution
In order to obtain the charged particle distribution (dN ch /dη) from dN/dη, either different bin-by-bin corrections can be used, or neutral particles can be marked as background ("fake") using RooUnfold's Fake() method. In this analysis, the following methods were used as the charged factor correction: 1. Bin-by-bin correction of the already unfolded dN/dη using the charged particle fraction N charged (η)/N all (η) from MC data; 2. Bin-by-bin correction of the raw EPD data via N charged (i Ring )/N all (i Ring ) from MC data; then unfolding of the EPD charged particle distribution. 7 3. Mark neutral particles as background and fill the response matrix as in the second method, except that the hits from neutral primaries are considered as "fake".
The three different methods can later be used to estimate the systematic uncertainty of the unfolding procedure itself.

Consistency check of the unfolding methods
Before unfolding the real data, a closure test was done to check whether the unfolding method can recover the "true" training data itself (MC "truth").
It was found that unfolding done on the input training MC sample reproduces well the input η distribution. When some noise (±1-10%) was added to the training sample, the resulting Figure 4: Consistency check of the three different methods to get dN ch /dη from MC EPD ring distribution. The difference is shown as unfolded dN ch /dη over MC "truth", the distributions divided bin-by-bin. Blue marker represents the first method (η-dependent charged factor correction), black shows the second method (EPD ring number dependent charged factor correction), and red represents the third method (marking neutral particles), relative to MC truth's dN ch /dη. The errorbars are only plotted for informative purposes: they were calculated using the ROOT's TH1 class' default square root of sum of squares of weights.
unfolded distribution was in agreement with the input distribution within < 3%. All in all, the unfolding itself was found to work well.
Furthermore, after applying the multiple counting correction and the three different methods of charged factor correction on the unfolded distribution 8 , the resulting distributions were compared to the original MC dataset's dN ch /dη. As it is visible in Fig. 4, the maximal relative deviation is up to 2% in certain bins for the first method and less than 0.1% for the other two methods.
It is worth noting that although the third method (marking neutral particles) shows here the most precise result, the systematic checks showed that it is the least reliable, in terms of most heavily depending on the MC input provided to the response matrix.
Given the result of the closure test, the unfolding and correction methods were considered adequately self-consistent.

Systematic checks
In the following section, the examined systematic uncertainty sources and their contribution to the results are discussed.

Dependence on input MC distribution
The bayesian iterative unfolding process, via its iterative nature, should mostly overcome differences in response matrix from real response that are not related to distortion effects, such as detector geometry [8]. However, as the exact response matrix cannot be determined even with precise MC simulations and the unfolding process itself is not perfect, some dependencies on the various parameters in the MC simulations can occur. Those are considered as systematic uncertainties of the measurement.

Tightening and shifting the input MC dN/dη
Firstly, the simulated sample's dN ch /dη was modified ("suppressed") using a Gaussian shape with width σ and mean η 0 . These suppression factors can be seen in Fig. 5a. This was done via a random selection based on Gaussian distribution while filling the response matrices.
Using this approach, all combinations could be analysed, that is, unfolding the i-th MC sample's EPD ring hit distribution via response from j-th MC sample. In case of i = j, the unfolding was as close to perfect as expected, discussed in Subsec. 2.4.
Unfolding results with the Gaussian width of σ 1 were not considered here as in this case there are almost no particles in the EPD range. Otherwise, there was less than a few percent Overall, in the analysis the effect of tightening the dN ch /dη of the training sample to σ = 2 and shifting it by ±3 units of pseudorapidity was investigated.

Broadening the input MC dN/dη
Similar to modification done in Sec. 3.1.1., here the tracks were modified with a factor of exp There was no suppression utilized for |η| > η max , with η max = 6. The resulting shape of the distributions can be seen in Fig. 6. While unfolding the data with these input MC distributions, a significant decrease at midrapidity values was observed. However, this occurred mostly outside the EPD's η region; the unfolding was considered acceptable down to σ broad ≈ 3.

Changing the charged fraction in the MC training dataset
The fraction of the charged particles in the MC input data was changed by ±15%. This was achieved by randomly rejecting either the neutral or the charged particles.

Changing the p T slope of the MC training dataset
The transverse momentum (p T ) distribution slope of the MC input data was changed by ±10% via randomly rejecting particles of small or large p T .

Centrality and z-vertex selection
It was investigated, by how much the unfolded distribution would change if either the z-vertex or the centrality selection are modified. For the former investigation, a ±5 cm calibration uncertainty in the z-vertex measurement of the real EPD data was employed; for the second one, ±5% calibration uncertainty was assumed in centrality determination of the real EPD data.

z-vertex choice
Due to the detector geometry, it is important to also take into account the interaction point's z-vertex position in the calculations, as the resulting pseudorapidity distribution should not depend on it.
The EPD data, as well as the responses, were collected in nine different z-vertex classes, equally distributed from −45 to +45 cm. Depending on which range was unfolded, the resulting distribution still may differ and has to be taken into account as systematic uncertainty.

Unfolding method choice
The most significant systematic uncertainty contribution was caused by the difference between the results achieved using different unfolding and correction methods (as listed in Subsection 2.3.). The first method was used as benchmark, from which the differences were calculated.

EPD related uncertainties
As previously stated, the EPD electronics were considered fully efficient (except some "dead areas" in the detector from e.g. glue and gaps, but these were assumed to be correctly handled in the simulation). The uncertainty from multi-MIP Landau fit was considered negligible compared to other systematic sources.
In conclusion, the systematic uncertainties coming from the detector system itself were considered negligible.
The different systematic uncertainty sources and their contribution with informative percentage values can be seen in Table 1.

Results
In this manuscript, charged particle pseudorapidity distributions with systematic uncertainties listed in Sec. 3. were obtained at two RHIC energies, in the EPD pseudorapidity range. The results at √ s N N = 19.6 and 27.0 GeV can be seen in Fig. 7 and 8, respectively. The caption #MIP ≤ 5 written on the plot indicates the number of convolution members in the multi-MIP Landau fit, as described in Sec. 1.1.

Comparison with the PHOBOS results
Another experiment of the RHIC complex was the PHOBOS experiment, which completed data taking in 2006. The PHOBOS was a large acceptance silicon detector, covering almost 2π in azimuth and |η| < 5.4 in pseudorapidity [9]. Compared to STAR's EPD, there are differences in both detector topology and granularity: the silicone pad detectors measure the total number of charged particles emitted in the collision, with modules mounted onto a centrally located octagonal frame (Octagon) covering |eta| ≤ 3.2, as well as three annular frames (Rings) on either side of the collision vertex, extending the coverage out to |η| ≤ 5.4 [12]. The PHOBOS also measured dN ch /dη at 19.6, 62.4, 130, 200 GeV energies [13]. Although in that paper a slightly different centrality binning was used (0-3%, 3-6% and 6-10% instead of 0-5% and 5-10%; the other centrality classes were the same), at 19.6 GeV the results can be compared.
In Figure 9, it is apparent that the two measurements show sizeable differences, depending on η: around up to a factor of two, increasing from small |η| towards forward/backward rapidities. The exact reasons behind this discrepancy are not yet known but the difference cannot be explained by the systematic uncertainties described in Sec. 3.

Discussion
In summary, based on EPD ring-by-ring distributions, charged particle pseudorapidity measurements at √ s NN = 19.6 and 27.0 GeV were performed with detailed systematic investigations regarding simulation data, calibration data, and unfolding methods. The results at √ s N N = 19.6 GeV show significant difference compared to the results from PHOBOS. There are four components in this comparison: EPD spectrum measurement, Geant4 simulation, unfolding procedure from the STAR part, and the PHOBOS data itself. The method presented in this manuscript is to be extended to other √ s N N values (as part of the BES-II program) and to fixed target data -mainly at energies where the QCD critical point is expected [14]. Refining this measurement method is also important for the search of the QCD critical point, in order to fine-tune the models used in these analyses.
Measuring pseudorapidity values of charged particles is important due to the possibility of estimating the initial energy density of the quark-gluon plasma created in the collisions, based on them [15,16]. Furthermore, the forward and backward rapidity measurements can provide information about the nuclear-matter effects as well [17].