First measurement of the sixth order cumulant of√
net-proton multiplicity distributions in sNN = 200 GeV
Au+Au collisions at the STAR experimentToshihiro Nonaka

February 2018

First measurement of the sixth order cumulant of√
net-proton multiplicity distributions in sNN = 200 GeV
Au+Au collisions at the STAR experimentToshihiro Nonaka
Doctral Program in Physics

Submitted to the Graduate School ofPure and Applied Sciencesin Partial FulFillment of the Requirements
for the Degree of Doctor of Philosophy in
Scienceat theUniversity of Tsukuba

Abstract
The existence of quarks and gluons has been explained by the standard model. In our universe,they are confined in a hadron and never observed solely due to the asymptotic freedom and theconfinement explained by the Quantum Chromo Dynamics (QCD). A few micro seconds after the Big
Bang, however, it is predicted that quarks and gluons are deconfined, which is called the Quark Gluon
Plasma (QGP). The understanding of QGP would possibly reveal the nature of the early universe.
Heavy ion colliding experiments have been carried out in order to form the QGP in the laboratory.
Based of the various results obtained in the last two decades, we have definitely succeeded to form the
QGP in the heavy ion colliding experiments. The hadron gas phase and the QGP phase are thus drawn
in the QCD phase diagram as a function of the temperature (T) and the baryon chemical potential(µB ). The next step is to understand the phase structure between those two phases. According to
the Lattice Gauge Theory calculation, the transition from QGP to hadronic system at small µB has
been considered as a smooth crossover, while lots of model calculations predict the first order phasetransition at large µB . If the first order phase transition exists at large µB , the QCD critical end
point should exist at the lower density end of the phase boundary. In Relativistic Heavy Ion Collider(RHIC), the Beam Energy Scan (BES) program is performed to search for the QCD phase structure.
By changing the colliding energy, various different T-µB regions would be explored in the QCD phase
diagram. Several recent results from the BES program the STAR experiments indicate possible signal
from the first order phase transition and the critical point, although they are not yet conclusive. Besidethese intensive searches for the phase structure at large µB , experimental inputs for understanding thestructure at small µB are still missing, where a smooth crossover is predicted by the Lattice Gauge
Theory calculation. The large number of observables are expected to be not very much sensitive tothe crossover transition due to the predicted continuity of the phase transition.
Fluctuation of conserved quantities is one of the powerful tools to study the QCD phase structure
especially to find the critical end point, where the correlation length could change or diverge. One of thepossible ways to test the prediction for a smooth crossover is to measure the higher order cumulantsof net-baryon or net-charge multiplicity distribution. Net-proton multiplicity distributions can bestudied as a reasonable proxy for net-baryon distributions. The STAR experiment has measured
the fourth order cumulant ratio (κσ 2 = C4 /C2 ) of net-proton multiplicity distribution in Au+Au√
collisions with its value of about ∼ 0.92 at sNN = 200 GeV, which is consistent with the model
prediction of a hadronic gas. Generally the higher order the cumulant, the more sensitive it is to thechange of the correlation length. Therefore the signature of the phase transition could be observed
via the fluctuation measurements of the higher order cumulants. Correction methods on various
backgrounds and experimental artifacts have not been established yet. It is pointed out that theprevious methods used in the STAR experiment up to the fourth order cumulant would not be corrected
for the non-binomial detector effects, and possible additional effect from the volume fluctuations fromparticipants has not been fully considered and corrected. Since the higher order cumulant consists ofthe combination of all the lower order cumulants, therefore the sixth order cumulants will be largely
affected by those effects compared to the lower order cumulants. For the measurement of the sixthorder cumulant, it is thus essential to further investigate the existing method and to develop new
methods for various corrections for the final experimental and detector effects as well as for the initialphysics effects like the participant fluctuations.
More efficient methods have been developed for the efficiency correction on cumulants based on thebinomial model of the experimental detector efficiency. The calculation cost has been reduced by factorof 1/100 compared to the conventional previous formulas. A new unfolding method has been developed,which is necessary to correct the non-binomial experimental detector effect. The effect of the volume
fluctuation on higher order cumulants has also been investigated with two methods, the Centrality
Bin Width Correction (CBWC) and the Volume Fluctuation Correction (VFC). These methods are
implemented to the analysis on the sixth order cumulant of the net-proton multiplicity distribution√
in Au+Au collisions at sNN = 200 GeV. The centrality dependence of C6 /C2 shows systematicallynegative values from mid-central to the most central collisions, which might be interpreted as a signalof a smooth crossover. Results are compared to the Poisson and the binomial baselines. Over thewide centrality region, C6 /C2 shows systematically smaller values from the Poisson baseline. The
binomial baseline can describe the results better than the Poisson baseline, however results in 30–40
and 10–20% centralities show even smaller values than the binomial expectation. These experimental
results are also compared to the UrQMD model. It is found that the UrQMD model is closed to the
Poisson baseline. The results show systematic suppression compared to the UrQMD model for all
centralities except for the most central collisions. More statistics would be necessary in order to draw
more definite conclusions about the possible signals from the crossover phase transition.

IIAcknowledgement
I would like to appreciate Prof. ShinIchi Esumi for his kind guidance to the physics of the heavy ion
colliding experiment. I really enjoyed the discussions on fluctuation analysis, which was very exciting
so that I forgot about the time. My fluctuation analysis all started with Dr. Hiroshi Masui. He
kindly guided me to the fluctuation analysis in the STAR experiment. Especially, the automation
technique for the efficiency correction on cumulants was so useful and impressive. Prof. Xiaofeng Luo,Prof. Nu Xu and Prof. Bedangas Mohanty also helped me fit into the fluctuation analysis in STAR.
They kindly taught me the analysis details on the higher order fluctuations. Discussions with themwere so fruitful. I would like to express my special thanks to Dr. Masakiyo Kitazawa in Osaka University. His short-term class in Tsukuba University when I was a master course is one of the reason
for me to have an interested in the fluctuation analysis. His ideas on the efficiency correction were so
impressive and made the calculations quite easy and efficient. I would like to appreciate Dr. Takafumi Niida, Dr. Takahito Todoroki and Dr. Sanshiro Mizuno. They kindly supported me for our shift
taking, life in BNL and the analysis in the STAR environment.
My motivation on the heavy ion physics in Ph.D course was originally motivated by MRPC-TOF development during master course student. Prof. Tatsuya Chujo led me to the MRPC development andmade me realize the pleasure to contribute to the development of the fast timing detector. Prof. Motoi Inaba kindly provided me great supports on the electronics especially on pre-amplifiers for MRPC.
Prof. Yasuo Miake always surprised me by his shape questions, which became hints for further understanding of MRPC. Mr. Sumio Kato kindly managed the network environment.
I would like to thank Ritsuya Hosokawa for his kind support as Ph.D student in the same grade. I
was always surprised at his heavy drinking. I would like to appreciate Naoto Tanaka and Hiroshi Nakagomi. Fishing with them is a good memory. I would like to thank R. Aoyama and T. Sugiura. We
work together on the STAR experiment. Short stay with them in BNL for the shift taking was quite
hard. Technical and physics discussions with them were so exciting.
Many thanks to maser and doctoral students in our high energy nuclear collision group (utkhii) in
Tsukuba group, especially to H. Yokoyama, K. Gunji, T. Nakajima, S. Horiuchi, D. Watanabe, K.
Ohsima, H. Ozaki, K. Kihara, T. Kobayashi, W. Sato, M. Hirano, J. Lee, I. Sakatani, T. Shioya, H.
Yamamoto, K. Ito, Y. Fukuda, B. C. Kim, R. Koyama, S. Kudo, K. Matsunaga, Y. Rebaza, K. Sato,T. Ichisawa, H. Kato, D. Kawana, Y. Kawamura, K. Nakagawa and T. Suzuki. My reseach life and thecoffee time with them were quite comfortable. I would like to appreciate S. Sakai and Y. Watanabe
for the physics discussions and holding banquets.
I would like to give a special thanks to Hitomi Koyama, Ikuo Koyama, Kyoko Koyama, Ryota Koyama,
Yuta Koyama and Tora-hei Koyama for their continuous supporting and encouragement on my research life. They warmly welcomed me to their house, trip and fishing.
Finally, I would like to express my gratitude to my family, Yasuyuki, Kazuko and Takafumi, for theirunderstanding, scolding and encouraging for over two decades.
I wouldn’t be who I am now without all those who worked, played and lived with me. I will neverforget the feeling of gratitude and all the memories with them.
IContents
1 Introduction
1.1 Quantum Chromo Dynamics (QCD) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Asymptotic freedom and confinement . . . . . . . . . . . . . . . . . . . . . . .
1.1.2 QCD phase diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Search for the QCD phase structure . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Relativistic heavy ion collision . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Beam Energy Scan Program (BES) . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.3 Freeze-out line . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.4 Turn-off the QGP signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.5 First order phase transition . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.6 Critical end point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.7 Physics motivation – search for the crossover . . . . . . . . . . . . . . . . . . .
1.2.8 Technical motivation – investigate the analysis methods . . . . . . . . . . . . .
1.3 Moments and cumulants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 Additivity of cumulants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.3 Factorial moments and factorial cumulants . . . . . . . . . . . . . . . . . . . .
1.3.4 Cumulants of conserved quantities . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.5 Statistical baseline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2
2
2
3
3
4
4
5
5
6
7
8
10
11
11
12
13
13
14

2 Experiment
2.1 Relativistic Heavy Ion Collider (RHIC) . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 The STAR detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Time Projection Chamber (TPC) . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Time Of Flight (TOF) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Vertex Position Detector (VPD) . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Run by run QA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Event and track selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3 Proton identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17
17
17
18
20
22
27
27
27
29

3 Detector Effect
3.1 Efficiency correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Single efficiency bin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 Many efficiency bins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.3 Importance of precise efficiency correction . . . . . . . . . . . . . . . . . . . . .
3.1.4 Statistical error estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31
31
31
32
34
36

II

3.2

Unfolding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Binomial model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Embedding simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Single particle efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Efficiency distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37
37
44
46
46
47

4 Volume Fluctuation
4.1 Centrality Bin Width Correction (CBWC) . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Centrality Bin Width Correction (CBWC) . . . . . . . . . . . . . . . . . . . . .
4.1.2 Auto-correlation effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Volume Fluctuation Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Toy model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 CBWC vs VFC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Centrality resolution effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 UrQMD model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51
51
51
52
52
53
55
55
55
56

5 Results and Discussions
5.1 Systematic uncertainties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Up to 4th order cumulants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 6th order fluctuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4 Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5 Comparison with theoretical calculations . . . . . . . . . . . . . . . . . . . . . . . . . .

60
60
62
63
67
67

6 Conclusions

70

A Centrality definition

71

B Tips on the efficiency correction
B.1 Stirling number of the first kind . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.2 Factorial moments with the binomial model . . . . . . . . . . . . . . . . . . . . . . . .

74
74
74

C Net-proton distribution
C.1 Shape for the sixth order cumulant . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.2 Net-proton distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

76
76
77

D Robustness to the time dependence

80

3.3

III

List of Figures
1.1(Left) Coupling constant α as a function of momentum transfer [1]. (Right) Static
potential between quarks and antiquarks calculated by lattice QCD as a function of
distance [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Example of QCD phase diagram [3]. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 An image of the time evolution of the relativistic heavy ion collision [4]. . . . . . . . .
1.4 Extracted chemical freeze-out temperature Tch as a function of baryon chemical potential µB , compared with the theoretical predictions (left), and the kinetic freeze-out
temperature Tkin as a function of the average transverse radial flow velocity hβi (right)
measured at the STAR experiment in BES energies. [5]. . . . . . . . . . . . . . . . .
1.5 (Left) The charged hadron yields normalized by number of binary collisions as a functionof the number of participant nucleons at 3.0 < pT < 3.5 GeV/c and 4.0 < pT <4.5 GeV/c in BES energies [6]. (Right) The slope parameter of the directed flow at
midrapidity dv1 /dy|y=0 for proton, antiproton and netproton as a function of beamenergy [7]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.6 Published results of the ratio of higher order cumulants up to fourth order as a functionof beam energy for net-charge (left), net-proton (middle) and net-kaon (right) multiplicity distributions [8–10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.7 Recent STAR preliminary results on κσ 2 (left) and Sσ (middle) for net-proton multiplicity distributions measured at 0.4 < pT < 2.0 (GeV/c) as a function of beamenergy [11]. The theoretical prediction is also shown in the right hand side panel [12].
1.8 The energy density computed by lattice QCD as a function of the temperature [13]. .
B
1.9 (Left) The ratio of the sixth to second order baryon number susceptibility (χB
6 /χ2 )
as a function of temperature for vanishing and finite density µq /T . The green shadedBarea represents the chiral crossover region. (Right) The expected values of χB
4 /χ2 andB
. . . . . . . . . . . . . . . . . . . . .
χB
6 /χ2 with different freeze-out conditions [14].
1.10 An example of skewness and kurtosis [15]. . . . . . . . . . . . . . . . . . . . . . . . .
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9

An aerial photo of RHIC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Overview of the STAR detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A schematic of the STAR TPC [16] . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A schematic on the anode pad plane with one full sector [16]. . . . . . . . . . . . . .
The energy loss as a function of momentum with the Bethe-Bloch expectations. The
magnetic field is set to be 0.25 T [16]. . . . . . . . . . . . . . . . . . . . . . . . . . . .
√1/β as a function of momentum, measured in Run 9 at sNN = 200 GeV in p+p
collisions [17]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A two-side view on the MRPC detector [18]. . . . . . . . . . . . . . . . . . . . . . . .
A schematic on TOF trays. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The front view of one assembly (left) and a photograph of two assemblies (right) [19].
IV

3
4
5

6

7

8

9
9

10
12
18
19
21
21
22
23
23
24
25

2.10 The front view of one assembly (left) and a photograph of two assemblies (right) [19].
VPD
2.11 Correlation of the primary vertex position along the beam pipe measured by VPD Zvtx√
TPC
and measured by TPC Zvtx , in p+p collisions at sNN = 510 GeV recorded in 2012√
(left), and in Au+Au collisions at sNN = 200 GeV recorded in 2009 (right) [19].
VPD − Z TPC from which the Z
The insets show the difference between them ∆Z = Zvtx
vtx
vtx
resolution of VPD can be extracted. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.12 Averaged z-vertex position, radial vertex position, number of primary tracks, distanceof closest approach (DCA), transverse momentum, and pseudo-rapidity as a function of
run index in Run11. Mean values in the same trigger condition are shown in red solidlines, and ±3σ lines are shown green. . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.13 (a) Radial vertex position, (b) vertex position along the beam pipe, (c) the differenceof vertex positions along the beam pipe between TPC and VPD, (d) number of hit
points in TPC used for the track reconstruction (e) number of hit points in TPC used
for calculating dE/dx. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.14 (a) Number of primary tracks associated with TOF, and (b) the number of primary
tracks associated with TOF requiring β > 0.1 in the Run10 central trigger. (c) and (d)are for the Run11 minimum bias trigger. . . . . . . . . . . . . . . . . . . . . . . . . .
2.15 (Left) dE/dx measured by TPC as a function of momentum multiplied by charge.
(Right) m2 measured by TOF as a function of momentum multiplied by charge [20].
3.1
3.2

3.3

3.4

3.5

3.6

3.7

3.8

3.9

Comparison of number of necessary terms in correction formulas in two methods [21].
Deviation of the efficiency corrected values of cumulants using averaged efficiency, assuming (a) Gauss distribution and (b) distribution that has 5% smaller cumulants than
Poisson distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Efficiency corrected C6 and C6 /C2 and their statistical errors estimated by bootstrap,
as a function of 50 independent trials. Red dotted lines represent the averaged value
for 50 trials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Flowchart for our unfolding method by using simple toy model with Poisson distributions. (A) is the experimental true, (B) is the simulation true, (C) is the experimental
measured and (D) is the simulation measured distribution. RMrev is the reversed response matrix. (C.F. rec) and (C.F. gen) are correction functions in reconstructed andgenerated coordinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Cumulants up to fourth order as a function of iteration. Red lines represent the true
cumulants, black lines are incremental unfolding, and blue lines are conventional unfolding. Results of independent 30 trials are superimposed. . . . . . . . . . . . . . . .
Flowchart for our unfolding method by using critical a distribution. (A) is the experimental true, (B) is the simulation true, (C) is the experimental measured and (D) is thesimulation measured distribution. RMrev is the reversed response matrix. (C.F. rec)and (C.F. gen) are correction functions in reconstructed and generated coordinates. .
Two-dimensional distributions for protons and antiprotons at 1st , 20th and 80th iteration of incremental unfolding, and 200th iteration of conventional unfolding Initial
distributions for experiment and simulation are also shown. . . . . . . . . . . . . . . .
Net-proton distributions at 1st , 20th and 80th iteration of incremental unfolding, and200th iteration of conventional unfolding Initial distributions for experiment and simulation are also shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Cumulants up to fourth order as a function of iteration. Red lines represent the true
cumulants, black lines are incremental unfolding, and blue lines are conventional unfolding. Results of independent 30 trials are superimposed. . . . . . . . . . . . . . . .
V

25

26

27

28

29
30
35

36

37

39

40

41

42

42

43

3.10 Results of toy model calculation with 100 bootstrap samplings with 100k events. Cumulants up to fourth order are shown as a function of 100 independent trials. Red lines
represent the averaged value over trials. . . . . . . . . . . . . . . . . . . . . . . . . . .
3.11 An image for the implementation of efficiency bins in unfolding. Np(pbar), low(high) andnp(pbar), low(high) are generated and reconstructed number of (anti)protons at low(high)pT regions. Fp(pbar), low(high) is the fraction at low and high pT regions with respect tocombined pT . εp(pbar), low(high) is the efficiency in a corresponding bin. . . . . . . . . .
3.12 (Top) Unfolded distribution for protons and antiprotons at 0–5%, 5–10%, 30–40% and70–80% centralities. (Bottom) Cumulants up to fourth order as a function of centrality.
Blue squares are results of efficiency correction, and red circles are efficiency correction.
Any volume fluctuation correction is not applied. . . . . . . . . . . . . . . . . . . . .
3.13 TPC tracking efficiency (black), TOF matching efficiency (blue) and combined efficiency
at 0.8 < pT < 2.0 (GeV/c)(red) for protons as a function of pT for each centrality. . .
3.14 pT integrated efficiency for protons (open) and antiprotons (filled) at 0.4 < pT < 0.8(GeV/c) (red) and 0.8 < pT < 2.0 (GeV/c) (blue) as a function of refmult3. Dottedlines are polynomial fittings. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.15 TPC tracking efficiency for protons at 0–5% centrality as a function of azimuthal angles
and rapidity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.16 Hypergeometric and beta-binomial distributions with N = 10 and ε = 0.6. Distributions with different α are shown in different colors. Black points represent a binomial
distribution. Bottom panels show the ratio to binomial. . . . . . . . . . . . . . . . . .
3.17 Number of reconstructed protons are plotted with different number of embedded protons
from 7 to 13 in 0–5% centrality at 1.0 < pT < 2.0 (GeV/c), which are fitted by binomial(red), hypergeometric (blue) and beta-binomial (green) distributions. Fitting resultsfor χ2 /NDF and α are also shown in each panel. . . . . . . . . . . . . . . . . . . . . .
3.18 χ2 /NDF extracted from Fig. 3.17 as a function of number of embedded protons forbinomial (red), hypergeometric (blue) and beta-binomial (green) fittings. . . . . . . .
3.19 The non-binomial parameter α extracted from the beta-binomial fitting to the embedding data as a function of the number of input (anti)protons. Data points are fitted by
y = a/x + b to extrapolate to the large number of (anti)protons. . . . . . . . . . . . .
4.1
4.2
4.3

4.4

4.5

κσ 2 (C4 /C2 ) as a function of centrality with different rapidity coverage for centrality
definition using UrQMD model [22]. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
κσ 2 (C4 /C2 ) of net-proton distribution as a function of centrality with different centrality definition using UrQMD model [22]. . . . . . . . . . . . . . . . . . . . . . . . .
Upper panel shows the net-proton distributions at 0–5, 20–30 and 60–70% centralities.
Npart fluctuates according to the Glauber model for red lines, while Npart is fixed inblack lines. Lower panels show the ratio between two net-proton distributions. . . . .
Cumulants up to fourth order as a function of centrality. For black circles Npart fluctuates which includes participant fluctuation given by the Glauber model For red squares
Npart is fixed, which are true value of cumulant in this model. Black solid lines are
results of VFC, and blue crosses are results of CBWC. . . . . . . . . . . . . . . . . . .
(Left) Multiplicity distributions with 100% (red), 50% (blue) and 20% (green) centrality
resolution generated by Glauber and two component model. Distributions with 50%
and 20% centrality resolution were defined by random sampling from the distribution
with 100% resolution. (Right) Results up to fourth order cumulant as a function ofcentrality without CBWC nor CBWC, with CBWC and with VFC from left to right.
Dotted lines represent the true cumulants. . . . . . . . . . . . . . . . . . . . . . . . .
VI

44

45

45
46

47
47

48

49
50

50

52
53

55

56

57

4.6

4.7
4.8

(a) The refmult3 distribution and (b) the Npart distribution, where the colors represent
the different centralities. (c) and (d) show the two-dimensional histogram for Np andNpbar at Npart = 240 in 10-20 and 0-10% centralities. . . . . . . . . . . . . . . . . . .
Cumulants up to sixth order as a function of centrality for CBWC-N (true), CBWC-M
and VFC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C3 /C2 , C4 /C2 and C6 /C2 as a function of centrality for CBWC-N (true), CBWC-M
and VFC. The ratios of CBWC-M and VFC to CBWC-N are shown in bottom panels.

Centrality dependence of C6 /C2 in Run11 without any volume fluctuation corrections.
Results from the efficiency correction are shown in red circles, and results from theunfolding with the beta-binomial model are shown in green crosses. . . . . . . . . . .
5.2 Cumulants up to fourth order as a function of centrality using the experimental data.
Results without CBWC nor VFC, with VFC and with CBWC. In the former two cases,
results with different centrality bin width of 8bin, 16bin and 32bin (corresponding to
2.5%, 5.0% and 10% centrality step) are shown. . . . . . . . . . . . . . . . . . . . . .
5.3 The centrality dependence of C3 /C2 normalized by the value of Skellam (top) and C4 /C2√(bottom) at sNN = 7.7, 19.6, 62.4 and 200 GeV. Results using CBWC are shown inblue crosses, and VFC results are shown in green triangles. The Skellam baselines areshown in dotted lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4 The STAR preliminary results on the beam energy dependence of κσ 2 = C4 /C2 and√Sσ/Skellam = (C3 /C2 )/Skellam using CBWC [23], where VFC results at sNN = 7.7,
19.6 and 62.4 and 200 GeV are superimposed in green circles, orange diamonds andpurple stars for 0-5, 5-10 and 70-80% centralities, respectively. . . . . . . . . . . . . .
5.5 Centrality dependence of C6 /C2 for Run10 and Run11. Statistical errors are shown inbars, while systematic uncertainties are shown in brackets. The expanded plot is shownin right hand side from 20% centrality. . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6 C6 /C2 as a function of centrality, where results from Run10 and Run11 are merged.
The expanded plot is shown in right hand side from 20% centrality. Skellam and thebinomial baselines are also shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.7 C6 /C2 as a function of centrality with Skellam baseline and the results of UrQMD
model. CBWC is applied in UrQMD model as is done in experimental data. . . . . .
5.8 C6 /C2 as a function of centrality using CBWC or VFC for both experimental data and
UrQMD model. Results with CBWC are shown in purple and the results with VFC
are in green. Results of UrQMD model are shown in band. . . . . . . . . . . . . . . .
5.9 Significance with respect to zero as a function of centrality for CBWC and VFC. The
expanded plot is shown in right hand side panel. . . . . . . . . . . . . . . . . . . . . .
B
5.10 χB
as a function of temperature [24]. The red
6 /χ2 calculated with the lattice QCD√
shaded area it the region of interest at sNN = 200 GeV. . . . . . . . . . . . . . . . .
√5.11 Final results on C6 /C2 of net-proton distribution at sNN = 200 GeV, where thecorresponding results calculated by the lattice QCD on the ratio of sixth to second
order susceptibility of the net-baryon distribution [24] are superimposed in red shadows.

58
59
59

5.1

A.1 Luminosity dependence of averaged refmult3 in −2 < V z < 0 cm. Blue markers areuncorrected and fitted by first order polynomial function (red lines) which is used forthe correction. Corrected results are shown in blue markers. . . . . . . . . . . . . . . .
A.2 (Left) Tail of the refmult3 distribution fitted by error function shown in a red line. The
blue dotted line represents the parameter c in Eq. (A.1). (Right) Maximum value of
refmult3 distribution as a function of Vz (black) and the corrected distribution (red). .
VII

61

62

63

64

65

65
66

66
67
68

69

71

72

1
A.3 (Left) Refmult3 distribution (black points) fitted by MC distribution generated by
Glauber and two-component model (red points). (Right) Ratio of MC to refmult3
distribution. The solid blue line represents the fitting for trigger inefficiency correction.
C.1 A simple sketch to stack n boxes for the n-th order cumulant. . . . . . . . . . . . . .
C.2 Gaussian distributions with the mean value of 0 and the sigma of 5, in which anothertwo symmetric Gaussians around N = ±15 are added (left) or subtracted (right). Redlines represent the Gaussian fittings. The ratio to the fitting are shown in the bottom
panels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.3 (Top) Net-proton distributions in 0-5, 30-40 and 70-80% centralities. Blue squares are
the measured distribution. Red circles are unfolded distribution without any volume
fluctuation corrections. Black lines represent Skellam distributions, and dotted green
lines are Gaussian distributions. (Bottom) Ratios of the unfolded distribution to the
Skellam and Gaussian distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.1 Relative deviations with respect to the true value of cumulants up to the sixth order asa function of the fraction of bad events in which the efficiency drops with ∆ε. . . . .

73
77

78

79

80

Chapter 1

Introduction
1.1

Quantum Chromo Dynamics (QCD)

Quantum Chromo Dynamics (QCD) is the theory of strong interactions between quarks and gluons.
Quarks and gluons have three kinds of color charge, red, green and blue, and for baryons and mesons
the color charge of white is required due to the combinations of constituent quarks. In this section, I
introduce the important characteristics of QCD.

1.1.1

Asymptotic freedom and confinement

In electrodynamics, the coulomb interaction between two particles is given bye1 e2 ~rF~ (~r) = k 2 , ~r = ~x − ~y , r = |~r|,
r r

(1.1)

where e denotes the charge of the particle and two particles are put in ~x and ~y . Equation (1.1)
indicates that the coulomb force becomes repulsive when the charges for two particles are identical,
and vice vera. It also indicates that the coulomb force becomes strong with small r = |~r| and becomes
weak with large r. In the case of QCD, however, the static potential Vs is given by
4 αs
Vs = − ×+ k × r,3r
12π
αs (µ2R ) ≈,(33 − 2nq )log(µ2R /Λ2QCD )

(1.2)
(1.3)

where αs is called running coupling constant with µR and ΛQCD being the QCD scale which are
determined experimentally. This αs becomes small in small r and becomes large in large r, whichis called asymptotic freedom. Asymptotic freedom was discovered by deep inelastic scattering (DIS)
experiments colliding a high energy electron to a proton, which gave an indirect evidence of quarks as
well. Equation (1.3) was shown analytically from the lowest term of perturbative expansion in highenergy scale. In low energy scale (long distance), αs becomes large, which could be the reason for theconfinement of quarks and gluons. However, higher order terms become more important in low energy
scale, thus, the confinement of quarks cannot be shown analytically. The confinement of quarks andgluons is shown numerically by lattice QCD [2]. Lattice QCD is the non-perturbative approach of thefirst principles. Infinite volume is approximated with finite space of lattice, and the QCD theory issolved on the lattice numerically.

2

CHAPTER 1. INTRODUCTION

3

Figure 1.1: (Left) Coupling constant α as a function of momentum transfer [1]. (Right) Static potential
between quarks and antiquarks calculated by lattice QCD as a function of distance [2].

1.1.2

QCD phase diagram

The state of the QCD matter can be characterized by the two-dimensional plane of the temperature(T ) and the baryon chemical potential (µB ), which is called the QCD phase diagram. Figure 1.2
shows one example of the QCD phase diagram. In high temperature and high density region, quarks
and gluons are released from the confinement, then the strongly interacting matter called quark gluon
plasma (QGP) is created. In high density and low temperature region, color superconductors are
predicted, which is also expected to exist in neutron stars. We note that the boundaries shown in
Fig. 1.2 are still predictions and there is no experimental evidence. Theoretically, it has been shown by
the lattice QCD that the phase transition between the hadronic phase and the QGP phase is smoothcrossover at vanishing baryon chemical potential [25, 26], while in nonzero µB , the first order phasetransition is predicted by various models. These predictions make it easy for us to expect the QCD
critical point [3, 27, 28], where the first order phase transition end.

1.2

Search for the QCD phase structure

For the last several decades, relativistic heavy ion collisions have been carried out to discover andstudy the QGP. Considering various experimental results, it may be no doubt that the QGP has beensuccessfully created experimentally. Our next step is to shed light on the QCD phase structure. So
far, however, we find only three things from the experiments as below:• The QCD phase diagram contains the QGP phase as well as the hadron gas phase.
• The chemical freeze-out line can be drawn in the phase diagram, which is very close to the chiralcrossover region estimated by the lattice QCD.
• The kinetic freeze-out line can be drawn at the low temperature side of the chemical freeze-out
line.

CHAPTER 1. INTRODUCTION

4

Figure 1.2: Example of QCD phase diagram [3].
Therefore, the QCD phase diagram in Fig. 1.2 is nothing more than an image expected from few
experimental inputs and so many theoretical predictions. Most importantly, experimental inputs on
the location and the order of the phase transition line are still missing. In this section, we introduce
some experimental attempts and results to reveal the QCD phase structure following the dynamics on
the relativistic heavy ion collisions.

1.2.1

Relativistic heavy ion collision

A few micro seconds after the Big Bang, our universe would be the QGP. The high energy heavy ion
colliding experiment is the unique way to create the QGP and study the initial stage of our universein the laboratory. Fig. 1.3 shows an image of the time evolution of the relativistic heavy ion collision.
In the initial stage with two accelerated heavy ions, color glass condensate (CGC) could be realizedin each heavy ion. At the collision of two heavy ions, hard scattering occurs, and the thermalizationis achieved. The QGP can be formed if the colliding energy is substantially large. At this time, thehydrodynamical expansion starts. With decreasing the temperature and density of the system, thehadronization of QGP starts. After hadronization ceases, the chemical freeze-out occurs in which
the yields of various hadron species are fixed, followed by the kinetic freeze-out, after when no more
momentum transfer occurs.

1.2.2

Beam Energy Scan Program (BES)√
The beam energy scan program has been carried out at sNN = 7.7, 11.5, 14.5, 19.6, 27, 39, 62.4 and
200 GeV in 2010, 2011 and 2014 at RHIC in order to search for the QCD phase structure. Heavy
ion collisions with high beam energy allows nucleons to slip through each other with dropping large
energy into the system, which leads to high temperature and low density of the created medium.
Once we decrease the beam energy, however, the stopping effect of nucleons becomes dominant. Thenthe created system will be low temperature with high density. Therefore, the structure of the QCD
phase diagram can be scanned and investigated by changing the colliding beam energy. From the next

CHAPTER 1. INTRODUCTION

5

Figure 1.3: An image of the time evolution of the relativistic heavy ion collision [4].
subsection, we introduce some results from BES and discuss the expected scenarios on the QCD phasestructure.

1.2.3

Freeze-out line

Two kinds of freeze-out are considered in heavy ion collisions: the chemical freeze-out and the kinetic
freeze-out. The chemical freeze-out is the time when inelastic collisions cease and particle yields
are fixed. Thus, the chemical freeze-out parameters are extracted by fitting the particle yields withthe statistical thermal model. After the chemical freeze-out, elastic interactions among particles willcontinue until the average interparticle distance becomes large enough, then momenta of particles
becomes unchanged, which is called the kinetic freeze-out. The kinetic freeze-out parameters can beobtained by the simultaneous fit to the pT spectra for different particles using the blast wave model.
Figure 1.4 shows the extracted chemical freeze-out temperature Tch as a function of baryon chemical
potential µB (left), and the kinetic freeze-out temperature Tkin as a function of the average transverse
radial flow velocity hβi (right) measured at the STAR experiment in BES energies. The value of Tkin issmaller than Tch . This is consistent with the picture of the freeze-out, the chemical freeze-out occurs
followed by the kinetic freeze-out. The separation between Tkin and Tch increases with increasing
the beam energy, which might be due to the increase of hadronic interactions between the chemical
and kinetic freeze-out at higher beam energies. The centrality dependence of µB is observed, whichbecomes strong in low beam energies, while Tch is similar for all centralities. Importantly, the lattice
QCD suggests the phase transition temperature around Tc ≈ 150 − 160 MeV, which is close to thechemical freeze-out temperature. It may indicate that the phase transition occurs right before thechemical freeze-out.

1.2.4

Turn-off the QGP signal

When we search for the existence or locations of the first order phase transition and the QCD criticalpoint by using certain observables, we have to implicitly require that the QGP is formed at thebeam energy. If not, we should not find any signals for the phase transition. In order to study
onset of charged particle suppression relative to the binary collision scaling, the normalized charged
hadron yields per binary collision have been measured at two pT bins of 3.0 < pT < 3.5 GeV/c and
4.0 < pT < 4.5 GeV/c in BES energies [6], as shown in the left hand side panel in Fig. 1.5. If the
QGP is formed, high pT partons lose their energy by passing though the volume of the QGP, andthe hadron yields decrease with respect to low pT . This is called the jet-quenching effect. Thus, thesuppression of high pT hadron yields has been considered to be an evidence of the QGP being formed.
On the contrary, however, an absence of this effect does not rule out the possibility of the QGP,CHAPTER 1. INTRODUCTION

6

Figure 1.4: Extracted chemical freeze-out temperature Tch as a function of baryon chemical potentialµB , compared with the theoretical predictions (left), and the kinetic freeze-out temperature Tkin as afunction of the average transverse radial flow velocity hβi (right) measured at the STAR experimentin BES energies. [5].
since there are other effects leading to the enhancement of hadron yields such as the Cronin effect,
radial flow and particle coalescence. At 7.7 GeV and 11.5 GeV the yields increases monotonically withincreasing hNpart i. This would not prove that there is no jet-quenching effect due to the absence of the
QGP, but simply that the enhancement effects grow rapidly than jet-quenching effect. The fact, thatthe yields become flat in central collisions 14.5 GeV, might be understood that the enhancement andsuppression effects are comparable. In higher energies the yields show suppression in central collisions,and 200 GeV data shows the monotonic decrease with increasing hNpart i, which would be because thejet-quenching effect is stronger than other enhancement effects. We note again that the observed√
enhancement of yields at the beam energy sNN ≤ 14.5 GeV does not exclude the scenario that the
QGP is still created in such low beam energy region. But considering that other measurements using
low pT probes also show some turn-off signal of the QGP in similar beam energy regions, the scenario√of the vanishing QGP phase at the beam energy of sNN ≤ 14.5 GeV may not be a bad estimate as
all [29, 30].

1.2.5First order phase transition

One possible observable to probe the first order phase transition is the directed flow v1 , which is thefirst harmonic coefficient of the Fourier expansion of the azimuthal distribution. It has been predicted√by hydrodynamic models that the dv1 /dy shows a minimum near midrapidity around sNN = 4 GeV
for net-baryons [31]. The right hand side panel in Fig. 1.5 shows the slope parameter of the directed
flow at midrapidity dv1 /dy|y=0 for proton, antiproton and netproton as a function of beam energy [7].
The antiproton slope shows monotonic decrease with decreasing the beam energy, whereas the proton
slope shows non-monotonic behaviour with the minimum value between 11.5 GeV and 19.6 GeV.
This cannot be explained by the UrQMD model. Furthermore, net-proton slope exhibits sign change
twice between 7.7 –11.5 GeV and between 19.6 –39 GeV. These results on the proton and net-proton
slope are qualitatively consistent with hydrodynamics predictions, which might be a possible signal

CHAPTER 1. INTRODUCTION

7

for the first order phase transition. But we need to take care that the beam energy in which thepredictions suggest the minimum value for the slope parameter is different from the observation.
Further theoretical understandings are necessary.

Figure 1.5: (Left) The charged hadron yields normalized by number of binary collisions as a functionof the number of participant nucleons at 3.0 < pT < 3.5 GeV/c and 4.0 < pT < 4.5 GeV/c in BES
energies [6]. (Right) The slope parameter of the directed flow at midrapidity dv1 /dy|y=0 for proton,
antiproton and netproton as a function of beam energy [7].

1.2.6

Critical end point

One of the most important results for the last ten years is the higher order fluctuation of conserved
quantities, which provides us some possible signals for the location of the QCD critical point. It
has been suggested that the measurements of higher order cumulants or moments (see Sec. 1.3) forthe event-by-event net-charge, net-baryon and net-strangeness multiplicity distributions as a functionof beam energy would show a non-monotonic behaviour with respect to the statistical baseline, if
the trajectory of the system created in heavy ion collisions pass through near the critical point [32].
Figure 1.6 shows the published results of σ 2 /M = C2 /C1 , Sσ = C3 /C2 and κσ 2 = C4 /C2 as afunction of beam energy for net-charge, net-proton and net-kaon multiplicity distributions at theSTAR experiment [8–10]. Net-proton and net-kaon are proxies for net-baryon and net-strangeness,
respectively. Results of net-charge and net-kaon distributions show mostly monotonic behaviour with√respect to the statistical baseline, while a non-monotonic kink is observed at sNN = 19.6 GeV for
net-proton. More detailed studies on the net-proton fluctuations has been carried out [11], in which
the non-monotonic behaviour of κσ 2 = C4 /C2 is observed as shown in Fig. 1.7 by extending the pT
coverage from 0.4 < pT < 0.8 GeV/c to 0.4 < pT < 2.0 GeV/c. Since this non-monotonic behaviour

CHAPTER 1. INTRODUCTION

8is qualitatively consistent with the model calculation [12], it is considered as a possible signal fromthe QCD critical point. However, there are still large statistical errors in the low beam energy region.
More statistics is required in order to derive more definite physics messages, which would be a main
goal of the Beam Energy Scan Phase II starting from 2019.

Figure 1.6: Published results of the ratio of higher order cumulants up to fourth order as a functionof beam energy for net-charge (left), net-proton (middle) and net-kaon (right) multiplicity distributions [8–10].

1.2.7

Physics motivation – search for the crossover

The phase transition at vanishing baryon chemical potential can be calculated by using the lattice
QCD, which has shown that the the phase transition is a smooth crossover without a clear boundary [26]. Figure. 1.8 shows the energy density computed by lattice QCD as a function of temperature [13]. It can be found that the energy density increases rapidly but smoothly starting fromT ≈ 150 MeV. However, there is no experimental evidence for this smooth crossover, because allobservables would not show any discontinuity for the crossover. Theoretically, the model calculation
predicts the sixth order cumulant for net-charge and net-baryon multiplicity distributions would be√come negative in sN N ≥ 60 GeV. The left hand side panel in Fig. 1.9 shows the ratio of the sixthB
to second order baryon number susceptibility (χB
6 /χ2 ) as a function of the temperature for vanishingBand finite density µq /T = 0, 0.14 and 0.44. For all the cases, we see the negative value of χB
6 /χ2 near
the chiral crossover region 0.9 < T /Tpc ≈ 1.0 and more than unity in T /Tpc < 0.9. The right handB
B
B
side panel in Fig. 1.9 shows the expected values of χB
4 /χ2 and χ6 /χ2 with different freeze-out conditions [14]. As discussed in SubSec. 1.2.3, the crossover region might be close enough to the chemical
freeze-out line. Thus, the physics motivation of this study is to find the signal for the crossover phasetransition by the measurement of C6 /C2 for the net-proton multiplicity distribution as a function of√centrality at sNN = 200 GeV. When considering the negative value of C6 /C2 at the phase transition
region, having more than unity for C6 /C2 would provide us another scenario of the crossover region
located far away from the chemical freeze-out line.

CHAPTER 1. INTRODUCTION

9

Figure 1.7: Recent STAR preliminary results on κσ 2 (left) and Sσ (middle) for net-proton multiplicity
distributions measured at 0.4 < pT < 2.0 (GeV/c) as a function of beam energy [11]. The theoreticalprediction is also shown in the right hand side panel [12].

Figure 1.8: The energy density computed by lattice QCD as a function of the temperature [13].

CHAPTER 1. INTRODUCTION

10

B
Figure 1.9: (Left) The ratio of the sixth to second order baryon number susceptibility (χB
6 /χ2 ) as afunction of temperature for vanishing and finite density µq /T . The green shaded area represents theB
BBchiral crossover region. (Right) The expected values of χB
4 /χ2 and χ6 /χ2 with different freeze-out
conditions [14].

1.2.8Technical motivation – investigate the analysis methods

The experimental measurement of higher order fluctuations of conserved quantities has a relatively
short history. The first publication up to the fourth order cumulants was carried out in year 2010
from the STAR experiment without the efficiency correction (see Sec. 3.1). The second publicationfrom STAR was done in 2014 with the efficiency correction to correct the finite detector efficiency.
The Centrality Bin Width Correction (see Sec. 4.1) was applied in order to correct for a possible
volume fluctuation in both publications. There were no other experimental groups until the HADES
experiment presented the results up to the fourth order cumulants for the net-proton distribution in
QM2017. In HADES preliminary results, their analysis methods were quite different from those inSTAR. They tried unfolding technique to correct the detector effects including efficiency. In addition,
they used another method (see Sec. 4.2) to correct the volume fluctuation instead of the Centrality Bin
Width Correction. In this way, developments for the analysis technique on higher order fluctuations
have been carried out by various people/groups/experiments. In other words, however, there arestill no established ways to correct these effects. Furthermore, it is expected that the sixth orderfluctuation are much more affected by these effects than the fourth order fluctuation, since the higher
order cumulant consists of all the lower order cumulants or moments as will be seen in Eqs. (1.8)–(1.20). In order to accomplish the physics motivation mentioned above, the technical development isthus essential. Therefore, the existing and newly developed analysis methods have been investigatedby applying to several models and the experimental data in order to extract more definite physics
messages and to look for some scenarios including the signal from the phase transition.

CHAPTER 1. INTRODUCTION

11

1.3

Moments and cumulants

1.3.1

Definitions

Moments and cumulants are mathematical measures which represent the shape of the distribution.
rth order moment µr are defined as:
drµr = hmr i = r G(θ) θ=0 ,dθXmθG(θ) =e P (m) = hemθ i,

(1.4)
(1.5)

m

P
where P (n) denotes a probability distribution function satisfying n P (n) = 1, G(θ) is a moment
generating function, and a bracket represents the average over many events. Similarly, rth order
cumulant Cr is given by
Cr = hmr ic =drK(θ) θ=0 ,dθr

K(θ) = lnG(θ),

(1.6)
(1.7)

where K(θ) denotes a cumulant generating function. As is seen in Eqs. (1.4)–(1.7), moments andcumulants are connected each other via generating functions. Thus moments are expressed in termsof cumulants asµ 1 = C1 ,

(1.8)

µ2 =

C2 + C12 ,

(1.9)

µ3 =

C3 + 3C2 C1 + C13 ,C4 + 4C3 C1 + 3C22 + 6C2 C12 + C14 ,C5 + 5C4 C1 + 10C3 C2 + 10C3 C12 + 15C22 C1 + 10C2 C13 + C15 ,C6 + 6C5 C1 + 15C4 C2 + 15C4 C12 + 10C32 + 60C3 C2 C1 + 20C3 C13+15C23 + 45C22 C12 + 15C2 C14 + C16 ,(1.10)

µ4 =
µ5 =
µ6 =

(1.11)(1.12)(1.13)and vice versa,
C1 = µ 1 ,(1.14)µ2 − µ21 ,(1.15)

C3 = µ3 − 3µ2 µ1 + 2µ1 ,(1.16)

C2 =C4 =µ4 − 4µ3 µ1 − 3µ22 + 12µ2 µ21 − 6µ41 ,(1.17)C5 =µ5 − 5µ4 µ1 − 10µ3 µ2 + 20µ3 µ1 + 30µ22 µ1 − 60µ2 µ31 + 24µ51 ,µ6 − 6µ5 µ1 − 15µ4 µ2 + 30µ4 µ21 − 10µ23 + 120µ3 µ2 µ1 − 120µ3 µ31+30µ32 − 270µ22 µ21 + 360µ2 µ41 − 120µ61 .

(1.18)

C6 =(1.19)

Due to the complicated expression in higher order moments and cumulants, it is convenient to use thefollowing recursion formula to implement the calculation:

r−1 
X
r−1
Cr = µr −Cn µr−n ,n−1n=1(1.20)

CHAPTER 1. INTRODUCTION

12

where we can calculate rth order cumulant automatically by taking appropriate combinations for up
to r − 1th order moments.
Generally, moments shown in Eqs. (1.5) and (1.8)–(1.13) are called non-central moments. This is often
compared to central moments. The rth order central moment µ0r are defined as the parameters withrespect to the mean value, which is given byµ0r = hδmr i, δm = m − hmi,

(1.21)

where µ01 = 0 by definition. Central moments are used for intuitive understanding of the shape of the
distribution. The second order central moments is called standard deviation σ, the third order moment
is the skewness S, and the four order moment is kurtosis κ. The standard deviation represents thewidth of the distribution. The skewness and the kurtosis represent the asymmetry and the sharpnessof the distribution, respectively. The positive of negative value of the skewness and the kurtosis is
illustrated in Fig. 1.10. The positive skewness gives the asymmetric shape leaned to the right handside with respect to the mean position of the distribution, whereas the shape becomes left hand side
leaned shape for the positive skewness. For the positive kurtosis, the center of the distribution becomes
sharper and the tail becomes wider. The negative kurtosis provides the box-like shape.

Figure 1.10: An example of skewness and kurtosis [15].

1.3.2

Additivity of cumulants

Although cumulants and moments are compatible measures each other, cumulants are more easy tohandle due to their additivity [15]. Let us consider two independent probability distribution functions
P1 (m1 ) and P2 (m2 ), where m1 and m2 are stochastic variables. Probability distribution function ofm = m1 + m2 are then given byX
P (m) =
δm,m1 +m2 P1 (m1 )P2 (m2 ).
(1.22)
m1 ,m2

Moment and cumulant generating functions for P (m) are calculated as follows:
X
X
X
G(θ) =emθ P (m) =emθ
δm,m1 +m2 P1 (m1 )P2 (m2 )
m

=

X
m1

m

e

m1 θ

P1 (m1 )Xm1 ,m2m2 θe

P2 (m2 ) = G1 (θ)G2 (θ),(1.23)

m2

K(θ) = lnG(θ) = K1 (θ) + K2 (θ).

(1.24)

CHAPTER 1. INTRODUCTION

13

By taking r derivatives on Eq. (1.24), r-th order cumulant can be expressed as
hmr ic = hmr1 ic + hmr2 ic ,(1.25)

which indicates that any order of cumulant for the probability distribution functions for sum of two
independent stochastic variables are given by the sum of cumulants for probability distribution function
for each variable. On the other hand, moments are obtained by
hmi = hm1 i + hm2 i,
2(1.26)

hm21 i + hm22 i + 2hm1 ihm2 i,
hm31 i + hm32 i + 3hm21 ihm2 i + 3hm1 ihm22 i,

hm i =
3

hm i =(1.27)(1.28)

where correlation terms of moments between two probability distribution functions for each stochastic
variables appear. This is the reason why cumulants are preferred to moments for various usages. Notethat this ”additivity” of cumulants is only valid in the case of independent stochastic variables.

1.3.3

Factorial moments and factorial cumulants

For the efficiency correction discussed in Secs. 3.1.1 and 3.1.2, let us introduce factorial moments andfactorial cumulants that can also characterize probability distribution functions [15]. The nth orderfactorial moment are defined asdn
hmn if = hm(m − 1)...(m − n + 1)i = n Gf (s) s=1 ,(1.29)dsX
Gf (s) =sm P (m) = G(ln s),(1.30)mwhere Gf (s) denotes the factorial moment generating function. Similarly, the nth order factorial
cumulant is given bydn
Kf (s),dsn
Kf (s) = lnGf (s) = K(ln s),

hmn ifc =(1.31)(1.32)

where Kf (s) denotes a factorial cumulant generating function. For the convenience, the expression for
the factorial moments of net-particle distribution can be introduced by:
Fab =

∞ X
∞
X

P (m, m̄)

m=a m̄=bm!
m̄!
,(m − a)! (m̄ − b)!

(1.33)

where m and m̄ are number of positively and negatively charged particles.

1.3.4

Cumulants of conserved quantities

The cumulants of conserved quantities in quantum statistical mechanics are considered. The grandpartition function for conserved quantity N is defined as
Z = tr[e−(H−µN )/T ],(1.34)

where H is Hamiltonian, V is volume, T is temperature, and µ is chemical potential. Cumulants of
N are given by derivatives of the grand potential Ω = −T ln Z with respect to µ/T :∂ n (−ω/T )V≡ χn V,∂(µ/T )n
Ω = ωVhN n ic =(1.35)(1.36)

CHAPTER 1. INTRODUCTION

14

where the grand potential per unit volume is defined as ω, and χn is defined as susceptibility. Experimentally, the cumulants hN ic can be calculated by looking at the event-by-event multiplicity
distribution for conserved quantities, e.g., net-charge, net-baryon and net-strangeness distributions.
Theoretically, susceptibilities can be computed by LQCD calculation. By taking ratio of different
order of cumulants, the volume V can be canceled, thus, experimental and theoretical results can bedirectly compared.

1.3.5

Statistical baseline

Binomial distribution
The binomial distribution is defined by the number of success in a sequence of N independent trials
with the success probability of p, which is given by
Bp,N (m) =N!
pm (1 − p)N −m .
m!(N − m)!

(1.37)

This is equal to the case that we measure the number of particles m with probability p in certain
volume with respect to the total number of N particles in the whole volume. The generating functions
are thus obtained as
GB (θ) =Xemθ Bp,N (m) =N!
(eθ p)m (1 − p)N −m
m!(N
−
m)!
m

X

m

= (1 − p + eθ p)N ,(1.38)

θ(1.39)KB (θ) = N ln(1 − p + e p).

By taking derivatives on the cumulant generating function, cumulants up to sixth order are given by
hmn ic = ξn N,(1.40)

withξ1 = p,(1.41)

ξ2 = p(1 − p),(1.42)

ξ3 = p(1 − p)(1 − 2p),(1.43)
2

ξ4 = p(1 − p)(1 − 6p + 6p ),(1.44)
2

ξ5 = p(1 − p)(1 − 2p)(1 − 12p + 12p ),23(1.45)
4

ξ6 = p(1 − p)(1 − 30p + 150p − 240p + 120p ).

(1.46)

Now one can consider the binomial baseline for the net-proton distribution by assuming protons andantiprotons follow the binomial distributions. The measured mean µ and the scaled variance σ 2 µ forthe distribution of proton and antiproton are written byµ = N p,ε=σ2= 1 − p,µ(1.47)

CHAPTER 1. INTRODUCTION

15

Cumulants up to fourth order can be expressed by using µ and ε:hmic = µ,(1.48)

2

hm ic = µε,(1.49)

3

hm ic = µε(2ε − 1),
4(1.50)

2

hm ic = µε(6ε − 6ε + 1),
52(1.51)

2

hm ic = µε (2ε − 1)(12ε − 12ε + 1),
6

4

3

2

hm ic = µε(120ε − 240ε + 150ε − 30ε + 1).

(1.52)(1.53)

Then the rth order cumulant of the net-proton distribution is given by
Crnet = Crp + (−1)r Crp̄ .

(1.54)In the experiment, measured C1 and C2 for the (anti)proton distribution allow us to calculate Eqs. (1.48)–(1.53), which provides the binomial baseline for the net-proton distribution as shown in Eq. (1.54).
Poisson distribution
By taking p → 0 limit of the binomial distribution, one obtains the Poisson distribution:
Pλ (m) =

λm λ
e ,
m!

(1.55)

where λ denotes the mean value of the distribution. The cumulant generating function is given byKλ (θ) = λ(eθ − 1).

(1.56)Taking derivatives on this equation provides cumulants:
hmn ic = λ (n ≥ 1).

(1.57)

It can be found that the all orders of cumulants are identical, and they are exactly equal to the meanvalue λ. Next, let us consider the difference of two Poisson distributions, which corresponds to the
net-proton distribution. The probability distribution function is written asX
Sλ1 ,λ2 (m) =δm,m1 −m2 Pλ1 (m1 )Pλ2 (m2 ),(1.58)
m1 ,m2

which is known as the Skellam distribution. Then the generating functions are given by
X
X
G(θ) =emθδm,m1 −m2 Pλ1 (m1 )Pλ2 (m2 )m=X

m1 ,m2em1 θPλ1 (m1 )m1

Xe−m2 θ P( λ2 )(m2 )m2= Gλ1 (θ)Gλ2 (−θ),K(θ) = Kλ1 (θ) + Kλ2 (−θ).

(1.59)(1.60)

By taking derivatives on the cumulant generating function, the rth order cumulant is given by
hmr i = λ1 + (−1)n λ2 ,(1.61)

CHAPTER 1. INTRODUCTION

16

which indicates that the odd order cumulant is equal to the difference of two means of two Poisson
distributions, whereas the even order cumulant is given by the sum of them. Therefore, the Skellam
baseline for cumulant ratios for the net-particle distribution are written as
R2,1 =

λ1 + λ2
C2=,C1
λ1 − λ2

R3,2 =

C3
λ1 − λ2=,
C2
λ1 + λ2

R6,2 = R4,2 = 1.

(1.62)

Experimentally, the Poisson distribution expresses the probability distribution to measure the numberof particle m in a substantially small volume. If one decreases the acceptance smaller and smaller tothe limit for measuring protons and antiprotons, the cumulants should become close to the values of
the Poisson distribution. In this case, there will be no chance to be sensitive to the fluctuation arisingfrom the possible critical point or phase transition. Keeping the acceptance as large as possible is thus
important in the fluctuation analysis, and comparing the measured fluctuations with respect to the
Skellam baseline would provide us some information on whether the observed signal is anything morethan the statistical fluctuation or not.

Chapter 2

Experiment
2.1

Relativistic Heavy Ion Collider (RHIC)

Search for the QGP started from the late 80’s in AGS and SPS accelerator in BNL and CERN by
colliding heavy ions to the fixed target. Although some signs of hot and dense matter were observed√at the SPS top energy of sNN = 17 GeV, it was not a definitive evidence. Relativistic Heavy
Ion Collider (RHIC) is the world’s first heavy ion collider starting from the year of 2000. Figure 2.1
shows an aerial photo of RHIC. RHIC consists of two concentric rings called ”blue ring” and ”yellow
ring” with circumferential length of 3.8 km, made up of 1,740 superconducting magnets. Two rings
intersect at six points, where different detectors were located at four of the intersection points, STAR,PHENIX, BRAHMS and PHOBOS. BRAHMS and PHOBOS completed their operations about ten
years ago, and PHENIX is being dismantled since last year. STAR is the only experiment running in
RHIC. At RHIC, various kinds of beam ion species can be accelerated, from the light nuclei such as√
deuteron and He3 to the heavy ions like Cu, Au and U. The highest beam energy is sNN = 510 GeV√
in p+p collisions and sNN = 200 GeV in Au+Au collisions.

2.2

The STAR detector

The Solenoid Tracker At RHIC (STAR) is the detector located at the 6 o’clock in RHIC. Figure 2.2
shows an overview of the STAR detectors. STAR has large and uniform acceptance in midrapidity,
which covers the pseudorapidity of −1.8 < η < 1.8 and full azimuthal angles, with excellent particle
identification capability. As is seen in the name of STAR, the main tracker is Time Projection Chamber(TPC) which enables us three-dimensional tracking of produced charged particles. Due the large
amount of particles produced in heavy ion collisions, lots of information are lost once the particle
tracking is profiled into two-dimensional plane. Thus, the three-dimensional tracking achieved by TPC
is one of the important features in the heavy ion colliding experiment. The Time Of Flight (TOF)
detectors are located outside the TPC. TOF measures the timing when particles pass through thedetector with less than 100 ps resolution, which enables us the particle identification at the intermediate
momentum region. The Barrel Electro-Magnetic Calorimeter (BEMC) are located radially further
outside of TOF, which covers |η| < 1.0 with full azimuthal angles. BEMC is used to trigger highmomentum photons and electrons. The Endcap Electro-Magnetic Calorimeter (EEMC) of the same
system with BEMC is located at the forward rapidity 1.0 < η < 2.0 with the full azimuthal angles
as well. The Beam Beam Counter (BBC) sits in 3.5 m away from the collision point covering 3.3 <|η| < 5.0. The BBC consists of two rings of hexagonal scintillator tiles. It is used to measure relative
luminosity for polarized p+p collisions. The Zero Degree Calorimeter (ZDC) is placed 18 m away from17

CHAPTER 2. EXPERIMENT

18

Figure 2.1: An aerial photo of RHIC.
the center of TPC, downstream of the dipole magnets, covering < 2.5 mrad in the polar angle. The
ZDC is the hadronic calorimeter. It measures the energy of the spectator nucleons. The coincidence
signals from detectors on either side of the interaction region are used for luminosity monitoring.
The Vertex Position Detector (VPD) is used to provide minimum-bias trigger in Au+Au collisions,
to measure the location of primary vertex along the beam pipe, and to determine the start time for
TOF.
From the next subsection, details for TPC, TOF and VPD, mainly used for fluctuation measurements,
will be explained.

2.2.1

Time Projection Chamber (TPC)

TPC is the primary tracking device in STAR, which performs tracking of charge particles to obtain
momenta and to identify particles by measuring the ionization energy loss in the gas volume. Momenta
are measured from 100 MeV/c up to 30 GeV/c. The schematic of TPC is shown in Fig. 2.3. The TPC
is 4.2 m long and 4 m in diameter, and it is located inside a solenoidal magnets with 0.5 T. The TPC
is filled with P10 gas (10% methane, 90% argon) regulated at 2 mbar above atmospheric pressure,
where the uniform electric field is applied ≈ 135 V/cm by Central Membrance, field-cage cylinders
and the readout end caps.
The readout system is based on Multi-Wire Proportional Chamber (MWPC). Once drifting electrons reach the anode wires, avalanche occurs with the amplification gain of 1000–3000. Then image
charge is induced by moving the large amount of positively charged ions. Anode wires are locatedwith 4 mm spacing, and the position resolution is thus limited by the spacing in this direction. The
direction of anode wires are set to be perpendicular to the tracks with the highest transverse momentum. Momentum resolution of charged particles gets worth with increasing pT since the radiusof the curvature of the track becomes large. As MWPC has the best position resolution in the wire
direction, better position resolution for higher pT charged particles can be achieved by setting theCHAPTER 2. EXPERIMENT

Figure 2.2: Overview of the STAR detector

19

CHAPTER 2. EXPERIMENT

20

anode wire direction roughly perpendicular to the straight tracks. Figure 2.4 shows a schematic on
the anode pad plane of one full sector. For the outer radius subsectors, readout pads are occupied
with no space between pad rows in order to optimize the dE/dx resolution. By collecting more signals
from the ionization electrons, statistics on the dE/dx measurement can be improved. For the inner
radius subsectors, smaller pad size is required due to the high track density. This is important for
charged particles with lower momentum passing across the pad rows far from perpendicular. Compromising with the limited number of readout electronics, separated pad rows are used for the inner
radius subsectors instead of the continuous pad coverage.
The x and y position of electron clusters are determined by the induced signals in a single pad row. The
signals are induced only on three pads for typical tracks. By assuming that the signals are distributedwith Gaussian, x position and its width with pad h2 centered at y = 0 is given by
 
σ2
w2
h3
x=
, σ2 =
,
(2.1)
ln
2w
h1
ln(h22 /h1 h3 )where h1 , h2 and h3 are amplitudes on three pads and w is the pad width. This algorithm gets worsein the case of the large crossing angles, when any three adjacent pads have comparable amplitudeof signals. Thus the weighted mean algorithm is also used as well. The z position of the track is
determined by measuring the time of secondary electrons reaching the anode wire to induce signals.
The z position is calculated by dividing the time by the drift velocity of electrons. This indicates thatthe reconstruction of z position strongly depends on the drift velocity, and it can change by the smallchange of the gas composition and the atmospheric pressure. In order to consider these effects, thedrift velocity is measured every few hours by artificial tracks created using the reference laser beams.
The main purposes of TPC is the vertex reconstruction, the measurement of momentum and theidentification of particles. The vertex resolution in the transverse plane is improved by the square rootof the number of tracks used in the vertex reconstruction. A resolution is achieved to be about 350 µm
with more than 1000 tracks. The momentum resolution is estimated by the embedding simulation,
where the simulated tracks from Geant are embedded in the real event and the tracking is performed.
Looking at the difference between true and reconstructed momentum of embedded tracks provides therelative momentum resolution. The best momentum resolution for pions is 2 % at pT ≈ 400 MeV/c,and 3 % for antiprotons at pT ≈ 800 MeV/c. The energy loss of charged particles in TPC gas volume
are calculated by summing up the induced charged on the readout pads. After considering the changeof the gas gain and the uncertainties arising from electronics, the dE/dx is extracted from the energy
loss measured at up to 45 pad rows. Figure 2.5 shows the energy loss as a function of momentum.
The magnetic field is set to be 0.25 T. It can be found that the pions and protons are separated up
to p ≈ 1 GeV/c.

2.2.2

Time Of Flight (TOF)

As is seen in Fig. 2.5, the particle identification between pion and proton with TPC is limited inp ≤ 1 GeV/c. The particle identification at the intermediate momentum region more than p ≈ 1 GeV/ccan be achieved by using additional information from Time Of Flight (TOF) detectors. TOF measures
the time when charged particles pass through with good resolution. As the distance L between thecollision point and the TOF is known, the mass is given by
 

t 2
2
2
m =p
−1 ,
(2.2)
L
where t is the time difference between the start and stop time. The stop time is given by TOF, andthe start time is measured by VPD (will be discussed in the next subsection). Figure 2.6 shows theCHAPTER 2. EXPERIMENT

Figure 2.3: A schematic of the STAR TPC [16]Figure 2.4: A schematic on the anode pad plane with one full sector [16].

21CHAPTER 2. EXPERIMENT

22

Figure 2.5: The energy loss as a function of momentum with the Bethe-Bloch expectations. The
magnetic field is set to be 0.25 T [16].
inverse of β = v/c as a function of the momentum provided by TPC. It can be found that pions,
kaons, protons and deuterons are clearly separated in this momentum region, for example, proton andpion separation below 3 GeV/c.
The TOF system in STAR is based on the Multi-gap Resistive Plate Chamber (MRPC). A schematicon the side view of MRPC is shown in Fig. 2.7. Six small gas gaps of 220 µm are created with highresistive glass plates and fishing lines. They are stacked between two electrodes on which high voltages±7 kV are applied. Then the high and uniform electric fields are formed in each small gas gap. Oncea charged particle pass through the MRPC, gas molecules are ionized to produce secondary electrons
and ions. Due to the high electric fields in a gas gap, an avalanche starts to grow up and move to theanode direction. This occurs at each gas gap but stops when an avalanche reach an adjacent glass
plate. Image charges are induced on the readout pad by fast moving of avalanche electrons. Signals
induced by secondary ions are negligible compared to those from electrons because they move slower
than electrons due to their large mass. Independent six avalanches are produced simultaneously withrespect to the incident of one charged particle, where the superposition of these independent signals
would lead an improved timing resolution over the single gap case. Full installations of 120 trays for thebarrel TOF system were carried out in the year 2010 with the acceptance of |η| < 0.9 and full azimuthal√
angles, and around 80 ps time resolution was achieved in Au+Au collisions at sNN = 200 GeV.

2.2.3

Vertex Position Detector (VPD)

VPD measures the arrival time of tens to hundreds of photons from π 0 , and provides the information
on event triggering [19]. The identical assemblies are located at both east and west sides with respect
to the collision point, which corresponds to the pseudorapidity region of 4.24 ≤ η ≤ 5.1. The primary

CHAPTER 2. EXPERIMENT

Figure 2.6: 1/β as a function of momentum, measured in Run 9 atcollisions [17].

23√sNN = 200 GeV in p+p

Figure 2.7: A two-side view on the MRPC detector [18].

CHAPTER 2. EXPERIMENT

24

Figure 2.8: A schematic on TOF trays.
vertex position along the beam pipe, Zvtx is given by
Zvtx = c(Teast − Twest )/2,

(2.3)

where the Teast and Twest are the times measured by each assembly, and c is the speed of light. By
using these two times, the start time for TOF can be determined as
Tstart = (Teast + Twest )/2 − L/c,

(2.4)

where L is the distance between either assembly and the collision point. The front view of one assembly
and a photograph of two assemblies are shown in Fig. 2.9. VPD consists of nineteen detectors madeup by a scintillator and a PMT. Such well-segmented detector with nineteen readout channels provide
the good time resolution by the central limit theorem, and thus leads to the better primary vertex
resolution, which is given by√
σstart = σ0 / N ,
(2.5)
√
√
σ(Zvtx ) = (c/ 2)σ0 / N ,(2.6)
where N is the number of readout channels and σ0 represents the time resolution of a single readoutdetector. Figure 2.10 shows the time resolution of VPD at each readout channel with different collisions
and beam energies. It can be found that the time resolution becomes better with increasing the beamenergy. This is because VPD is doing multiple-particle timing, thus the time resolution gets betterin higher beam energy due to larger multiplicity. By taking average of the times measured at all
readout channels, the start time resolution reaches around 20 – 30 ps. The actual time resolution
for single particle is the quadruple sum of the time resolution for VPD (start) and TOF (stop), thus
less than 30 ps start time resolution of VPD is good enough compared to the stop time resolution of
TOF about 80 ps. Figure 2.11 shows the correlation of the primary vertex position along the beam
VPD and measured by TPC Z TPC , in p+p collisions at √spipe measured by VPD Zvtx
NN = 510 GeV√ vtx
recorded in 2012 (left), and in Au+Au collisions at sNN = 200 GeV recorded in 2009 (right). The
VPD − Z TPC from which the Z
insets show the difference between them ∆Z = Zvtx
vtx resolution of VPD
vtx
can be extracted, and it can be found that the resolutions are σ(∆Z) = 2.3 cm and 0.9 cm for p+p
and Au+Au collisions, respectively.

CHAPTER 2. EXPERIMENT

25

Figure 2.9: The front view of one assembly (left) and a photograph of two assemblies (right) [19].

Figure 2.10: The front view of one assembly (left) and a photograph of two assemblies (right) [19].

CHAPTER 2. EXPERIMENT

26

VPD
Figure 2.11: Correlation of the primary vertex position along the beam pipe measured by VPD Zvtx√
TPC
and measured by TPC Zvtx , in p+p collisions at sNN = 510 GeV recorded in 2012 (left), and in√
Au+Au collisions at sNN = 200 GeV recorded in 2009 (right) [19]. The insets show the difference
VPD − Z TPC from which the Z
between them ∆Z = Zvtx
vtx resolution of VPD can be extracted.
vtx

CHAPTER 2. EXPERIMENT

2.3

27

Datasets

The measurements of higher order fluctuations are very statistical hungry. In order to reduce statistical√errors, the datasets collected in 2010 and 2011 at sN N = 200 GeV are both analyzed independently
and merged. Details for the datasets are summarized in Tab. 2.1.
Year
2010
2010
2011

Triggerminimum bias (10–80%)
central trigger (0–10%)
minimum bias (0–80%)

Events
200 M
160 M
510 M

Table 2.1: Datasets used in the analysis.

2.3.1

Run by run QA

Figure 2.12 shows the averaged z-vertex position, radial vertex position, number of primary tracks,
distance of closest approach (DCA), transverse momentum, and pseudo-rapidity as a function of run
index in Run11. Outlier runs outside ±3σ at each trigger condition are rejected from the analysis.

Figure 2.12: Averaged z-vertex position, radial vertex position, number of primary tracks, distanceof closest approach (DCA), transverse momentum, and pseudo-rapidity as a function of run index in
Run11. Mean values in the same trigger condition are shown in red solid lines, and ±3σ lines areshown green.

2.3.2

Event and track selection

Event and track quality cuts are shown in Tab. 2.2, and also shown in Fig. 2.13. In order to reject pile
up events, difference between z-vertex measured by TPC and VPD are limited within 3 cm. Residual
pile up evens are excluded from the correlation between the number of primary tracks associated with
TOF and number of primary tracks in |η| < 1.0 measured by TPC (RefMult) as shown in Fig. 2.14-(a)and (c). Events including lots of wired tracks with wrong β information are rejected by looking at
the correlation between the number of primary tracks associated with TOF requiring β > 0.1 and theCHAPTER 2. EXPERIMENT

28

RefMult, which is shown in Fig 2.14-(b) and (d). Tracks associated to the primary vertex are selectedwith DCA<1.0 cm. Primary tracks with good tracking quality are selected by requiring hit points forfitting the hit points and calculating the dE/dx.
|Vz |Vr|V pdVz − Vz |tofmatched
tofmatchedbeta
tofmatched
tofmatchedbeta
DCA
nHitsFit
nHitsFit/nFitPoss
nHitsDedx< 30 cm< 2 cm< 3 cm
> 0.50 × refmult − 30 (Run10)> 0.43 × refmult − 15 (Run10)> 0.50 × refmult − 13 (Run11)> 0.46 × refmult − 10 (Run11)< 1 cm> 20> 0.52>5

Table 2.2: Datasets used in analysis.

Figure 2.13: (a) Radial vertex position, (b) vertex position along the beam pipe, (c) the difference ofvertex positions along the beam pipe between TPC and VPD, (d) number of hit points in TPC used
for the track reconstruction (e) number of hit points in TPC used for calculating dE/dx.

CHAPTER 2. EXPERIMENT

29

Figure 2.14: (a) Number of primary tracks associated with TOF, and (b) the number of primary
tracks associated with TOF requiring β > 0.1 in the Run10 central trigger. (c) and (d) are for the
Run11 minimum bias trigger.

2.3.3

Proton identification

Left hand side plot in Fig. 2.15 shows dE/dx measured by TPC as a function of momentum multipliedby charge. Protons and antiprotons are identified by using dE/dx measured by TPC at 0.4 < pT < 0.8(GeV/c). At 0.8 < pT < 2.0 (GeV/c), however, it is difficult to identify protons only by using TPC
due to the contamination from pions and kaons. Thus, at 0.8 < pT < 2.0 (GeV/c), combined PID
with TOF is performed by requiring 0.6 < m2 < 1.2 (GeV/c2 ), which is shown in Fig. 2.15.

CHAPTER 2. EXPERIMENT

30

Figure 2.15: (Left) dE/dx measured by TPC as a function of momentum multiplied by charge. (Right)
m2 measured by TOF as a function of momentum multiplied by charge [20].

Chapter 3

Detector EffectExperimentally, some particles with finite detector efficiency, which distorts the shape of event-by-event
net-proton distributions. Conventional way to correct this effect is the analytical efficiency correction
with the assumption of the binomial response of efficiency [33–36]. However, it is known that theefficiency correction does not work when the binomial assumption breaks. Unfolding is necessary tocorrect such effect. In this chapter, we show the conventional formulas for efficiency correction andtheir implementation. But the calculation cost for the conventional formulas increases with number of
efficiency bins. In order to solve this problem, we show more efficient formulas that are robust to theincreases of efficiency bins. Then we move to the discussions on details for our unfolding approach.
Finally we try to investigate the non-binomial effect by using the existing embedding datasets.

3.1

Efficiency correction

3.1.1

Single efficiency bin

They are all derived based on the binomial model,
B(n, N ; ε) =N!
εn (1 − ε)N −n ,n!(N − n)!

(3.1)

where N and n denote the number of produced and measured particles, and ε represents the efficiency.
When the efficiency follows the binomial distribution, the relationship between true and measured
factorial moments can be expressed asfi,j= εi εj Fi,j ,(3.2)

with F and f being true and measured factorial moments of the net-particle distribution. Since therelationships between measured and true cumulants are not straightforward, the efficiency correction
is carried out using Eq. (3.2) as following steps:(a) Convert measured cumulants to measured factorial moments(b) Convert measured factorial moments to true factorial moments using Eq. (3.2)(c) Express the true cumulants in terms of true factorial moments
Final correction formulas can be found in Ref. [33] up to the sixth order cumulant, but it is cumbersome
to implement those explicit formulas in our program. In order to automate the efficiency correction,
31

CHAPTER 3. DETECTOR EFFECT

32

we decompose factorial moments into the correlation terms:
fab =

X
a

s(a, i)M

b
X

i

i=0

Fab =

X
a

j

s(b, j)M=

j=0s(a, i)N

i

i=0

b
X

a X
b
X

j

s(a, i)s(b, j) M i M ,(3.3)

i=0 j=0s(b, j)N

j


=

i

a X
b
X

j

s(a, i)s(b, j) N i N ,(3.4)

i=0 j=0

j=0

thus

a X
b
X

M M

j

a b

=ε ε

i=0 j=0

a X
b
X

j

N iN ,(3.5)

i=0 j=0

where s(n, k) denotes the Stirling number of the first kind (see B.1). The derivation of Eq. (3.5) can
be found in B.2 Since s(i, i) = 1, we can deduce the following recursive expressions:
X
X
b
j
b
j
s(a, i)s(b, j) M i M = εa εb N a N + εa εbs(a, i)s(b, j) N i N ,M aM +i,j≥0
i,j6=a,b
a

→ N N

b

M aM
=
εa ε bi,j≥0i,j6=a,b
b

+



X

s(a, i)s(b, j)i,j≥0i,j6=a,b

M iM
εa εb

j


− N N
,
i

j

(3.6)

0

with the initial condition hn0 n0 i = hN 0 N i = 1. In the case of fourth order cumulant, for example,jwhat one needs to do is to calculate all combinations for M i M with i, j ≥ 0, where i and j cannotjbe 4 simultaneously. Then one obtains N i N (a, b ≤ 4) by taking appropriate combinations for
M iM

3.1.2

j

recursively as shown in Eq. (3.6).

Many efficiency bins

Experimentally, efficiency depends on pT , rapidity, azimuthal angles, and particle species, which needs
to be implemented in the correction formulas. Number of factorial moments that we need to calculate
for mth order cumulant with M efficiency bins is given by
fm
Nm
=

m
X

r+M −1 Cr =m+M Cm − 1,

(3.7)

r=1

, which indicates that the calculation cost becomes drastically large with increasing the efficiency bins.
Therefore, if we apply efficiency correction precisely with many efficiency bins, the calculation will
not finish with limited CPU time. In order to solve this problem, we have developed more efficient
formulas by using the simple relationship between true and measured factorial cumulants as is seen
in factorial moments in Eq. (3.2). We show final correction formulas up to sixth order cumulant andfourth order mixed cumulants.

CHAPTER 3. DETECTOR EFFECT

33

Q c = hq(1,1) ic ,

(3.8)

2
ic + hq(2,1) ic − hq(2,2) ic ,Q2 c = hq(1,1)

(3.9)

3
ic + 3hq(1,1) q(2,1) ic − 3hq(1,1) q(2,2) ic + hq(3,1) ic − 3hq(3,2) ic + 2hq(3,3) ic ,Q3 c = hq(1,1)(3.10)

222
4ic
q(2,2) ic + 4hq(1,1) q(3,1) ic + 3hq(2,1)
q(2,1) ic − 6hq(1,1)
ic + 6hq(1,1)
Q4 c = hq(1,1)
2
ic − 12hq(1,1) q(3,2) ic + 8hq(1,1) q(3,3) ic − 6hq(2,1) q(2,2) ic+3hq(2,2)+hq(4,1) ic − 7hq(4,2) ic + 12hq(4,3) ic − 6hq(4,4) ic ,(3.11)

2
2
3
3
5
q(3,2) ic
q(3,1) ic − 30hq(1,1)
q(2,2) ic + 10hq(1,1)
q(2,1) ic − 10hq(1,1)
ic + 10hq(1,1)
Q5 c = hq(1,1)
222+20hq(1,1)
q(3,3) ic + 15hq(2,2)
q(1,1) ic + 15hq(2,1)
q(1,1) ic − 30hq(1,1) q(2,1) q(2,2) ic+5hq(1,1) q(4,1) ic − 35hq(1,1) q(4,2) ic + 60hq(1,1) q(4,3) ic − 30hq(1,1) q(4,4) ic+10hq(2,1) q(3,1) ic − 30hq(2,1) q(3,2) ic + 20hq(2,1) q(3,3) ic−10hq(2,2) q(3,1) ic + 30hq(2,2) q(3,2) ic − 20hq(2,2) q(3,3) ic+hq(5,1) ic − 15hq(5,2) ic + 50hq(5,3) ic − 60hq(5,4) ic + 24hq(5,5) ic ,(3.12)

6443
3
Q6 c = hq(1,1)
ic + 15hq(1,1)
q(2,1) ic − 15hq(1,1)
q(2,2) ic + 20hq(1,1)
q(3,1) ic − 60hq(1,1)
q(3,2) ic
3
2
2
2
2
2
+40hq(1,1)
q(3,3) ic − 90hq(1,1)
q(2,2) q(2,1) ic + 45hq(1,1)
q(2,1)
ic + 45hq(1,1)
q(2,2)
ic
3
3
2
2
2
2
+15hq(2,1)
ic − 15hq(2,2)
ic + 15hq(1,1)
q(4,1) ic − 105hq(1,1)
q(4,2) ic + 180hq(1,1)
q(4,3) ic − 90hq(1,1)
q(4,4) ic
22−45hq(2,1)
q(2,2) ic + 45hq(2,2)
q(2,1) ic + 60hq(1,1) q(2,1) q(3,1) ic − 180hq(1,1) q(2,1) q(3,2) ic+120hq(1,1) q(2,1) q(3,3) ic − 60hq(1,1) q(2,2) q(3,1) ic + 180hq(1,1) q(2,2) q(3,2) ic − 120hq(1,1) q(2,2) q(3,3) ic+6hq(1,1) q(5,1) ic − 90hq(1,1) q(5,2) ic + 300hq(1,1) q(5,3) ic − 360hq(1,1) q(5,4) ic + 144hq(1,1) q(5,5) ic+15hq(2,1) q(4,1) ic − 105hq(2,1) q(4,2) ic + 180hq(2,1) q(4,3) ic − 90hq(2,1) q(4,4) ic−15hq(2,2) q(4,1) ic + 105hq(2,2) q(4,2) ic − 180hq(2,2) q(4,3) ic + 90hq(2,2) q(4,4) ic
222+10hq(3,1)
ic − 60hq(3,1) q(3,2) ic + 40hq(3,1) q(3,3) ic + 90hq(3,2)
ic − 120hq(3,2) q(3,3) ic + 40hq(3,3)
ic+hq(6,1) ic − 31hq(6,2) ic + 180hq(6,3) ic − 390hq(6,4) ic + 360hq(6,5) ic − 120hq(6,6) ic ,(3.13)

CHAPTER 3. DETECTOR EFFECT

34

hQ(x) Q(y) ic = hq(1,0,1) q(0,1,1) ic + hq(1,1,1) ic − hq(1,1,2) ic ,(3.14)

2
q(0,1,1) ic + 2hq(1,0,1) q(1,1,1) ic − 2hq(1,0,1) q(1,1,2) ic + hq(0,1,1) q(2,0,1) ic − hq(0,1,1) q(2,0,2) ic
hQ2(x) Q(y) ic = hq(1,0,1)+hq(2,1,1) ic − 3hq(2,1,2) ic + hq(2,1,3) ic ,(3.15)

22
ic
q(0,1,1)hQ2(x) Q2(y) ic = hq(1,0,1)2
2
22q(2,0,2) ic
q(2,0,1) ic − hq(0,1,1)q(0,2,2) ic + hq(0,1,1)q(0,2,1) ic − hq(1,0,1)+hq(1,0,1)+4hq(1,0,1) q(0,1,1) q(1,1,1) ic − 4hq(1,0,1) q(0,1,1) q(1,1,2) ic+2hq(1,0,1) q(1,2,1) ic − 6hq(1,0,1) q(1,2,2) ic + 4hq(1,0,1) q(1,2,3) ic+2hq(0,1,1) q(2,1,1) ic − 6hq(0,1,1) q(2,1,2) ic + 4hq(0,1,1) q(2,1,3) ic
2
2ic
ic + 2hq(1,1,2)−4hq(1,1,1) q(1,1,2) ic + 2hq(1,1,1)+hq(2,0,1) q(0,2,1) ic − hq(2,0,1) q(0,2,2) ic − hq(2,0,2) q(0,2,1) ic + hq(2,0,2) q(0,2,2) ic+hq(2,2,1) ic − 7hq(2,2,2) ic + 12hq(2,2,3) ic − 6hq(2,2,4) ic ,(3.16)3hQ3(x) Q(y) ic = hq(1,0,1)
q(0,1,1) ic
2
2+3hq(1,0,1)
q(1,1,1) ic − 3hq(1,0,1)
q(1,1,2) ic + 3hq(2,0,1) q(1,0,1) q(0,1,1) ic − 3hq(2,0,2) q(1,0,1) q(0,1,1) ic+3hq(1,0,1) q(2,1,1) ic − 9hq(1,0,1) q(2,1,2) ic + 6hq(1,0,1) q(2,1,3) ic+3hq(2,0,1) q(1,1,1) ic − 3hq(2,0,1) q(1,1,2) ic − 3hq(2,0,2) q(1,1,1) ic + 3hq(2,0,2) q(1,1,2) ic+hq(3,0,1) q(0,1,1) ic − 3hq(3,0,2) q(0,1,1) ic + 2hq(3,0,3) q(0,1,1) ic+hq(3,1,1) ic − 7hq(3,1,2) ic + 12hq(3,1,3) ic − 6hq(3,1,4) ic ,(3.17)

where we used the symbolq(r,s)MX= q(ar /ps ) =(ari /psi )ni ,i=1
MXq(r,s,t) = q(xr ys /pt ) =(xri yis /pti )ni .

(3.18)(3.19)

i=1

This work is published in Physical Review C [21], where we also discussed the importance of precise
efficiency correction with many efficiency bins.
We performed the calculation for efficiency correction in case of many efficiency bins up to 4+4, withsimple toy model for both factorial moment method and new method. It was found that both methods
give exactly the same value numerically. Comparison of number of necessary terms in correction
formulas are compared between two methods as a function of number of efficiency bins in Fig. 3.1. It
can be found that the number of terms are reduced by a factor of 100 in new methods N6c comparedto factorial moment method N6c , and also it does not depend on the number of efficiency bins. CPUtime is also compared in Tab. 3.1, where we see that the CPU time is drastically reduced comparedto the conventional formulas.

3.1.3

Importance of precise efficiency correction

In previous section, we show more efficient formulas for efficiency correction. Then the question is,
however, whether we really need to apply precise efficiency correction with many bins. In order toCHAPTER 3. DETECTOR EFFECT

35

Nfm
6

106

c

Number of terms

N6
Nfm
4

105

c

N4
4

10

103

M
4
8
12
200

factorial moment
64.7s
17.3 × 102 s
14.1 × 103 s–

new method
30.8s
31.3s
32.3s
62.7s

102
10
0

5
10
15
20
Number of efficiency bins (M)Table 3.1: Comparison of CPU time to calculate
the sixth order cumulant between the conventionaland new methods [21].

Figure 3.1: Comparison of number of necessary
terms in correction formulas in two methods [21].
answer this question, we performed simple analytical calculation. The setup is following. We assume
one distribution that we want to measure cumulants, which consists of two internal independent
distributions. Both distributions have exactly the same value of cumulants Cn , then the cumulants
for whole distributions is are 2Cn . Also we assume two distributions have different efficiencies. We
apply efficiency correction by using the averaged efficiency for two distributions, then the deviation∆Kn can be expressed as:
∆K2 =
∆K3 =
∆K4 =



1 ∆ε 2(C2 − C1 ),(3.20)2 ε
3 ∆ε 2(C3 − 2C2 + C1 ),(3.21)2 ε
1 ∆ε 2
1 ∆ε 4(6C4 − 18C3 + 19C2 − 7C1 ) +(C4 − 6C3 + 11C2 − 6C1 ), (3.22)
2 ε
8 εEquations (3.20–3.22) are plotted in Fig. 3.2 as a function of ∆ε with assuming the averaged efficiency
to be ε̄ = 0.5. Figure 3.2 (a) assumes Gausssian distribution where more than second order cumulants
are zero. Figure 3.2 (b) assumes 5% deviation of cumulants from the usual Poisson distribution. It canbe found that the deviations from the true values become large with increasing the efficiency difference,and higher order cumulants have large deviation compared to lower orders. Interestingly, ∆Kn becomeszero by substituting Cn = C1 to Eqs. (3.20)–(3.22) assuming the Poisson distributions. Since wemeasure the deviation from Poisson as a signal from the QCD critical point or phase transition, it isimportant for us to apply precise efficiency correction with appropriate efficiency bins.

CHAPTER 3. DETECTOR EFFECT

36

2

(a)

0

C1=20, C2=25, C3= 0, C4= 0
ε=0.5
∆ K2
∆ K3
∆ K4

100

∆K

∆K

200

(b)

−2
C1=20, C2=19, C3=19, C4=19

−4

ε=0.5

0
−6
0

0.2

0.4

∆ε

0.6

0

0.2

0.4

0.6

∆ε

Figure 3.2: Deviation of the efficiency corrected values of cumulants using averaged efficiency, assuming(a) Gauss distribution and (b) distribution that has 5% smaller cumulants than Poisson distribution.

3.1.4

Statistical error estimation

Two methods are currently suggested for statistical error estimation for higher order cumulants. One
is the Delta theorem and the other one is bootstrap [35, 37]. Delta theorem is the analytical way toestimate the statistical errors for cumulants and moments, but it is cumbersome to derive analyticalformulas for higher order cumulants. In addition, those formulas are much complicated after efficiency
correction, which will make the implementation very difficult in the case of many efficiency bins. On
the other hand, bootstrap is very simple. One makes the new net-proton multiplicity distribution
by random sampling from the original distribution and calculate cumulants. This is repeated witharound 300 times, then the RMS of cumulants with 300 independent trials is taken as the statisticalerrors. Any en route calculation such as CBWC and the efficiency correction can be included inside thebootstrap. Therefore, bootstrap is used for the statistical error estimation for sixth order cumulant.
In a simple case of single efficiency bin, bootstrap is performed by random sampling from the two dimensional histogram between protons and antiprotons, which can be easily done by TH2D::GetRandomfunction in ROOT. In the case of many efficiency bins, however, further implementations are necessary. Suppose we estimate the statistical errors in case of n+n efficiency bins. We prepare the 2ndimensional histograms which contains the bin coordinate with its contents. After normalization forbin contents, bin coordinates can be randomly selected by generating random numbers from 0 to 1
for each sampling.
In order to justify the bootstrap, we performed simple toy model calculations. The setup for the
toy model is following. We generate four Poisson distributions with 1M events, two are for protons andthe others are for antiprotons, where the mean parameter for Poisson is 5, 4, 3, and 1, respectively.
Those distributions are randomly sampled with corresponding efficiencies, 80%, 60%, 70% and 90%.
The we apply efficiency correction to calculate C6 and C6 /C2 . Similarly, this procedures are also
performed with respect to a new distribution made by random sampling from the original distribution,
which is repeated with 300 times with 300 different (but sampled from the same original distribution)
distributions. Then RMS of C6 and C6 /C2 are taken as the statistical errors These procedures are
randomly and independently performed with 50 times, which is shown in Fig. 3.3. The number ofpoints touching the averaged value for 50 trials is 37 and 32 for C6 and C6 /C2 . They correspond to

CHAPTER 3. DETECTOR EFFECT

37

74% and 64%, which are reasonable values of 1σ for a Gauss distribution to justify bootstrap for C6
and C6 /C2 .

Figure 3.3: Efficiency corrected C6 and C6 /C2 and their statistical errors estimated by bootstrap, asa function of 50 independent trials. Red dotted lines represent the averaged value for 50 trials.

3.2

Unfolding

Efficiency correction formulas are derived based on the binomial response of efficiency. This binomial
model can be broken in the experiment due to the track merging or splitting that occurs in the detector
level. Reference [38] firstly demonstrates the effect of non-binomial efficiency. Once we apply efficiency
correction under the non-binomial situation, then the results will be artificially changed, and the effectbecome large for higher order cumulants. Such non-binomial detector effect can no longer corrected byanalytical approach – we need unfolding. Unfolding is the method to correct the detector effect, which
is provided in ROOT as RooUnfold package [39]. However, most of programs are already packaged in
the library, and users cannot touch the details for the implementations. In addition, it is limited tothree dimensional unfolding. Since we should be totally clear for analysis methods especially to searchfor unknown distributions predicted around the QCD critical point, we start to develop new unfolding
methods from scratch. In recent work, it was found that our unfolding method can easily correct
the volume fluctuation, and it can be also extended to the multidimensional case, which is strong
advantage compared to the existing unfolding method. Unfortunately discussions here are limited tobasic part of our unfolding method, I hope those advanced unfolding method will be presented inother literature. Correction for non-binomial effect on higher order cumulants is first attempt in theSTAR experiment. unfolding for sixth order cumulant is world’s first attempt. In this section, weshow the methodology of our unfolding method, then consider some non-binomial scenario expectedin the STAR experiment from the embedding data.

3.2.1

Methodology

Poisson test
We explain the methodology with a simple toy model by using Poisson distributions. Figure 3.4
shows the flowchart for our unfolding method. As a setup, we first generate two two dimensional
histograms for protons and antiprotons, one is for (A) experimental true distribution and the other isCHAPTER 3. DETECTOR EFFECT

38

for (B) simulation true distribution. Since the simulations are not complete and there should be somedeviation compared to experiment, mean parameter for Poisson is set to be different between (A) and(B) as µexp,p = 10, µexp,pbar = 7, µsim,p = 12 and µsim,pbar = 9. We apply the detector effect which
we call ”MC filter”, here we simply assume binomial efficiency εp = 0.9 and εpbar = 0.7, to (A) and(B) to make (C) experimental measured and (D) simulation measured distributions. Our goal is toreconstruct (A) by using (B), (C), (D) and MC filter. Our method consists of two kinds of unfolding
approach. Let us start from the ”incremental” unfolding approach.
1. Define the correction functions by the difference between experimental and simulation measured
distributions (B)–(D) as shown in (C.F. rec).
2. Considering that such difference between experiment and simulation in reconstructed coordinate should also appear in generated coordinate, convert the correction function (C.F. rec) to(C.F. gen) by using the reversed response matrix RMrev .
3. Apply smoothing to the correction function (C.F. gen).
4. Apply the correction function (C.F. gen) with scaling factor α to the simulation true distribution(C) to get (C’) which is slightly modified to close to the experimental true distribution (D).
5. Apply MC filter again to the modified simulation true distribution (C’) to generate the modified
simulation measured distribution (D’), which is again compared to the experimental measured
distribution (B).
6. Repeat 2–5 with substantial iterations.
Reversed response matrix RMrev in 2 can be constructed in each iteration when we apply MC filter
to the simulation true distribution. Smoothing in 3 is required to make the distribution smooth.
Scaling factor α in 4 is essential to avoid the negative bin content after applying the correction
function. In this approach, we apply MC filter iteratively to the simulation true distribution. Startingfrom the any distribution in simulation true, it should be incrementally modified to the experimental
true distribution by comparing with the experimental measured distribution in each iteration. Thus
we call this approach ”incremental” unfolding. Cumulants in each iteration are plotted in Fig. 3.5.
Red lines represent the true cumulants and black lines are incremental unfolding. Many lines aresuperimposed, which represents 30 independent trials for Fig. 3.4. It can be found that the cumulants
of incremental unfolding becomes close to true values with iterations, but it seems difficult for C3 andC4 to recover the input values. This would be because we apply MC filter in each iteration, whichmight lead to unstable results of higher order cumulants. In order to solve this problem, another
unfolding approach is implemented that we call ”conventional” unfolding. After additional iterations
for incremental unfolding, we switch to the conventional unfolding where we don’t apply MC filter.
We let RMinc,for and RMinc,rev being the forward and reversed response matrix defined in the last
iteration for incremental unfolding.
1. Apply RMinc,rev to the experimental measured distribution to obtain the simulation true distribution.
2. Refold the simulation true distribution by using RMinc,for to make a simulation measured distribution.
3. Define the correction functions by taking the difference from the experimental measured distribution.

CHAPTER 3. DETECTOR EFFECT

39

4. After conversion of correction functions from reconstructed to generated coordinate, apply it to
the simulation true distribution.
5. Repeat 2–5 until cumulants for net-distribution converge.
Since we don’t apply MC filter and don’t update the response matrix, this ”conventional” unfolding
get quickly converged as shown in blue lines in Fig. 3.5.

Figure 3.4: Flowchart for our unfolding method by using simple toy model with Poisson distributions.
(A) is the experimental true, (B) is the simulation true, (C) is the experimental measured and (D) is thesimulation measured distribution. RMrev is the reversed response matrix. (C.F. rec) and (C.F. gen)are correction functions in reconstructed and generated coordinates.

Critical shape test
Since our new unfolding approach were originally motivated by reconstruction of the unknown shapeof the conserved charge distributions predicted around the QCD critical point, we performed a criticalshape test to show the feasibility of our unfolding method. Figure 3.6 shows the initial distributionsin the critical shape test. Critical shape is assumed to have larger (or wider) Gauss distribution than
the Poisson distribution at smaller (or larger) Np values. Our goal is to reconstruct this critical shapestarting from the simple Poisson distribution for simulation true distribution. Methodology is thesame with performed in Poisson test. 100 and 200 iterations have been performed for incremental
and conventional unfolding respectively. Two-dimensional distributions are shown in Fig. 3.7 at theCHAPTER 3. DETECTOR EFFECT

40

Figure 3.5: Cumulants up to fourth order as a function of iteration. Red lines represent the true
cumulants, black lines are incremental unfolding, and blue lines are conventional unfolding. Resultsof independent 30 trials are superimposed.

CHAPTER 3. DETECTOR EFFECT

41

different stage of iterations as well as initial distributions for experiment and simulation. A simple
Poisson distribution is modified step by step, and finally the critical shape has been recovered. Let
us also look at the net-distribution in Fig. 3.8. We see the two-peak structure for the critical shape.
Starting from the simple Skellam distribution, two-peak structure has been recovered. Resulting
cumulants are shown as a function of number of iterations in Fig. 3.9. We see that our unfolding
approach works well to give true cumulants, and it can be found that the value of C4 has beensuccessfully modified from the value of Skellam to -2500.

Figure 3.6: Flowchart for our unfolding method by using critical a distribution. (A) is the experimental
true, (B) is the simulation true, (C) is the experimental measured and (D) is the simulation measured
distribution. RMrev is the reversed response matrix. (C.F. rec) and (C.F. gen) are correction functions
in reconstructed and generated coordinates.

Statistical error
Statistical errors for unfolding can be estimated by bootstrap. A simple toy model calculation has
been performed to show the feasibility of bootstrap on the unfolding. Figure. 3.10 shows the cumulants
up to fourth order as a function of 100 independent trials. At each data point, bootstrap is performedby random sampling from the same experimental true distribution with 100 times. In the next trial,
another experimental true distribution are independently generated to which we apply bootstrap
sampling. This procedure is repeated with respect to 100 independent experimental true distribution.
It can be found that the probability of data points touching the averaged value over 100 trials arearound 70%, which shows that bootstrap works well for unfolding.

CHAPTER 3. DETECTOR EFFECT

42

Figure 3.7: Two-dimensional distributions for protons and antiprotons at 1st , 20th and 80th iteration ofincremental unfolding, and 200th iteration of conventional unfolding Initial distributions for experimentand simulation are also shown.

Figure 3.8: Net-proton distributions at 1st , 20th and 80th iteration of incremental unfolding, and 200th
iteration of conventional unfolding Initial distributions for experiment and simulation are also shown.

CHAPTER 3. DETECTOR EFFECT

43

Figure 3.9: Cumulants up to fourth order as a function of iteration. Red lines represent the true
cumulants, black lines are incremental unfolding, and blue lines are conventional unfolding. Resultsof independent 30 trials are superimposed.

CHAPTER 3. DETECTOR EFFECT

44

Figure 3.10: Results of toy model calculation with 100 bootstrap samplings with 100k events. Cumulants up to fourth order are shown as a function of 100 independent trials. Red lines represent the
averaged value over trials.

3.2.2

Binomial model

Our main motivation on the unfolding is to correct the non-binomial efficiency. Before we start tostudy the effect of non-binomial efficiency, however, it is important to check whether the unfolding
with binomial model gives the consistent results with efficiency correction or not. In the efficiency
correction, pT dependent efficiency has been implemented. For unfolding, it can be simply included
inside the MC filter as shown in Fig. 3.11. We know the fractions of low and high pT regions Fp(pbar), low
and Fp(pbar), high , which are calculated by efficiency corrected C1 (mean) value separately for low andhigh pT regions. Based on the information in a two-dimensional histogram in generated coordinate,
we can divide the number of (anti)protons to low and high pT regions. We apply the MC filter
with the binomial efficiency in a corresponding bin. They are merged to be filled into the twodimensional histogram in the reconstructed coordinate. Figure 3.12 shows the results of unfolding at√
sN N = 200 GeV. Unfolded distributions are shown in top panels. This is the first attempt to lookat the corrected distributions. Bottom panels are cumulants up to fourth order for both efficiency
correction and unfolding with binomial model. It can be found that both results are consistent, whichshows the validity of our unfolding method. Note that the statistical errors for unfolding are slightly
larger than efficiency correction, which is due to the MC statistics. Both statistical errors should be
identical if we have substantial CPU time to increase the MC statistics.

CHAPTER 3. DETECTOR EFFECT

45

Figure 3.11: An image for the implementation of efficiency bins in unfolding. Np(pbar), low(high) andnp(pbar), low(high) are generated and reconstructed number of (anti)protons at low(high) pT regions.
Fp(pbar), low(high) is the fraction at low and high pT regions with respect to combined pT . εp(pbar), low(high)is the efficiency in a corresponding bin.

Figure 3.12: (Top) Unfolded distribution for protons and antiprotons at 0–5%, 5–10%, 30–40% and70–80% centralities. (Bottom) Cumulants up to fourth order as a function of centrality. Blue squares
are results of efficiency correction, and red circles are efficiency correction. Any volume fluctuation
correction is not applied.

CHAPTER 3. DETECTOR EFFECT

46

3.3

Embedding simulation

3.3.1

Single particle efficiency

Tracking efficiency of TPC can be estimated by embedding simulation. Embedding simulation is theMC approach, where we embed protons into real experimental data in detector level to see whether theembedded tracks can be reconstructed or not. By applying the same cuts in analysis to the embedded
tracks, single track efficiency can be estimated. TOF matching efficiency are estimated according to
the following equation by using the real experimental data:
εTOF =|nσ| < 2 && TOF matched.
|nσ| < 2(3.23)

Figure 3.13 shows the TPC tracking efficiency, TOF matching efficiency and combined efficiency at
0.8 < pT < 2.0 (GeV/c) for protons as a function of pT for each centrality. Then we obtain the pT

Figure 3.13: TPC tracking efficiency (black), TOF matching efficiency (blue) and combined efficiency
at 0.8 < pT < 2.0 (GeV/c)(red) for protons as a function of pT for each centrality.
integrated efficiency by using the following expression:R
ε(pT )f (pT )pT dpTR,
εint =
f (pT )pT dpT(3.24)

where f (pT ) are taken from the corrected pT spectra for protons and antiprotons [40]. Figure 3.14
shows the pT integrated efficiency for protons and antiprotons at 0.4 < pT < 0.8 (GeV/c) and 0.8 <pT < 2.0 (GeV/c) as a function of refmult3. Those data are fitted with polynomial function, which
gives the value of refmult3 bin by bin efficiency.
In the net-proton analysis up to fourth order cumulant at the STAR experiment, four efficiency
bins (low and high pT , protons and antiprotons) have been implemented in the efficiency correction.
It was found that TPC tracking efficiency depends on rapidity and azimuthal angles as well as pT due
to the fact that some TPC sectors can be in bad condition with electronic issues. Figure 3.15 shows
the azimuthal and rapidity dependence of TPC tracking efficiency for protons in the 0–5% centrality.
We see that the efficiency in sectors −π/3 < φ < 0 are lower than others, which is also implemented
in this study. Then the number of efficiency bins becomes eight including azimuthal dependence.

CHAPTER 3. DETECTOR EFFECT

47

Figure 3.14: pT integrated efficiency for protons(open) and antiprotons (filled) at 0.4 < pT < 0.8 Figure 3.15: TPC tracking efficiency for protons at(GeV/c) (red) and 0.8 < pT < 2.0 (GeV/c) (blue) 0–5% centrality as a function of azimuthal angles
as a function of refmult3. Dotted lines are polyno- and rapidity.
mial fittings.

3.3.2

Efficiency distribution

Non-binomial distributions
By replacing the MC filter from binomial model to non-binomial model, we can see some scenarios
how final results can be changed. We introduce two kinds of non-binomial distributions discussed inRef. [41]. One is the hypergeometric distribution and the other one is the beta-binomial distribution.
Former one has narrower width than binomial distribution, while the latter one has wider width.
Definitions for these two distributions are following. Let us consider to determine the number ofmeasured particles n from the true particles N with efficiency ε according to the hypergeometric
distribution. Suppose an urn containing white and black balls. Number of those balls are denoted by
Nw and Nb . We draw a ball from the urn, and count up n if it is white. This procedure is repeatedwith N times without throwing back to the urn. On the other hand, the beta-binomial distribution
can be implemented similar manner but the previous procedure is repeated by returning two balls toNb
the urn. Efficiency is determined by number of white and black balls as 1−εε = Nw . The deviationfrom binomial distribution can be implemented by another parameter α which determines the numberof while balls. Figure 3.16 shows examples for hypergeometric and beta-binomial distributions. We
assumed N = 10, ε = 0.6 and Nw = 2αN . Different α parameters are shown in different colors.
Binomial distribution is shown in black points. It can be found that smaller α gives narrower(wider)width for hypergeometric(beta-binomial) distributions, and it becomes close to binomial distribution
with large α. This reflects the fact that the replacement or returning of balls will be irrelevant in thecase of large number of white and black balls.
Modeling of the efficiency
We can model the response function extracted from the embedding datasets by fittings with binomial
and non-binomial distributions, which will allow us to expect how much the efficiency can be nonbinomial. Figure 3.17 shows the results of fittings to embedding datasets in 0–5% centrality at 1.0 <pT < 2.0 (GeV/c) by using the binomial, hypergeometric and beta-binomial distributions. Each panel
shows different number of embedded protons, and the number of reconstructed protons are plotted in

CHAPTER 3. DETECTOR EFFECT

48

Figure 3.16: Hypergeometric and beta-binomial distributions with N = 10 and ε = 0.6. Distributions
with different α are shown in different colors. Black points represent a binomial distribution. Bottom
panels show the ratio to binomial.

CHAPTER 3. DETECTOR EFFECT

49

black points. Lines in different colors represent the different fitting functions. Results of χ2 /NDF are
shown in Fig. 3.18 as a function of number of embedded protons for different fitting functions. It canbe found that the beta-binomial distribution is the best function to describe the embedding datasets
followed by binomial and hypergeometric distributions. Fig. 3.19 shows the extracted α parameter bybeta-binomial fittings as a function of the embedded (anti)protons, which is fitted by y = a/x + b to
extrapolate to the large number of embedded (anti)protons. These parameter a and b will be used asa part of systematic studies.

Figure 3.17: Number of reconstructed protons are plotted with different number of embedded protons
from 7 to 13 in 0–5% centrality at 1.0 < pT < 2.0 (GeV/c), which are fitted by binomial (red),
hypergeometric (blue) and beta-binomial (green) distributions. Fitting results for χ2 /NDF and α are
also shown in each panel.

CHAPTER 3. DETECTOR EFFECT

50

Figure 3.18: χ2 /NDF extracted from Fig. 3.17 as a function of number of embedded protons forbinomial (red), hypergeometric (blue) and beta-binomial (green) fittings.

Figure 3.19: The non-binomial parameter α extracted from the beta-binomial fitting to the embedding
data as a function of the number of input (anti)protons. Data points are fitted by y = a/x + b to
extrapolate to the large number of (anti)protons.

Chapter 4

Volume Fluctuation
Experimentally, the volume fluctuation effect needs to be addressed as backgrounds. Even in one
centrality, two nuclei collides with different overlap volume or with different number of participantnucleons (Npart ) event by event. Then the number of produced particles fluctuates due to the volume
fluctuation, which directly affects the values of cumulants. In order to suppress the effect of the volume
fluctuation, the STAR experiment has been using Centrality Bin Width Correction (CBWC). Recently,
on the other hand, P. Braun-Munzinger. et al newly suggested a volume fluctuation correction (VFC)
to eliminate the participant fluctuation with some model assumptions [42]. In this chapter, we discussboth CBWC and VFC. We compare both methods by using simple toy models, and show some
important features.

4.1

Centrality Bin Width Correction (CBWC)

In this section, we note some important points on CBWC related to the centrality resolution and theauto-correlation effect using the UrQMD model [22]. Based on these points we can determine the best
condition to define the centrality with volume fluctuation being substantially suppressed. Centrality
definition is summarized in App. A.

4.1.1

Centrality Bin Width Correction (CBWC)

Centrality Bin Width Correction (CBWC) is defined as:XCr =ωi Kri ,

(4.1)

i

ωi =

n
Pi ,
i ni

(4.2)

where i denotes the ith multiplicity bin, Kri is the rth order cumulant in ith multiplicity bin, and ni is
the number of events in ith multiplicity bin. It can be found that CBWC depends on the multiplicity
distribution used for the centrality definition. In other words, CBWC will be changed by the centrality
resolution. Figure 4.1 shows κσ 2 of net-proton distribution as a function of centrality with different
rapidity coverage used for centrality determination. The value of κσ 2 is artificially enhanced withnarrow rapidity coverage for the centrality definition (bad centrality resolution), while κσ 2 is suppressed to converge to certain value with large coverage for the centrality definition (good centrality
resolution). This indicates that the better centrality resolution is required for better correction with
CBWC, and thus we need to measure as large number of particles as possible for better centrality
resolution better in the experiment.
51

CHAPTER 4. VOLUME FLUCTUATION

52

Figure 4.1: κσ 2 (C4 /C2 ) as a function of centrality with different rapidity coverage for centrality
definition using UrQMD model [22].

4.1.2

Auto-correlation effectFigure 4.2 shows the centrality dependence of κσ 2 for net-proton distribution with BES energies with
two sets of centrality definition, one uses all charged particles in |η| < 2.0 where protons used for
cumulant analysis are included in the centrality determination, while the other uses charged kaonsand pions in |η| < 2.0. We can see that the results for former case are suppressed compared to thelatter case due to the auto-correlation effect, which becomes large in low beam energies, because thenumber of protons commonly used for cumulant analysis and the centrality determination becomeslarge due to the baryon stopping. Note that due to the resonance decay there are still some fractionsof auto-correlation effect even if we define the centrality with pions and kaons excluding protons.

4.2

Volume Fluctuation Correction

As discussed in Sec. 4.1, CBWC has been studied with UrQMD model [22], where the conclusion is the
CBWC results converge to certain value by increasing the centrality resolution, where we need to keep
in mind the fact that CBWC results will not converge to the true cumulants even in UrQMD model.
Recently new method for the volume fluctuation correction has been suggested in Ref. [42]. Although
it needs some model inputs for correction factors, it has nice and important features compared to
CBWC. In this section, we first derive the correction formulas. Then we perform a simple toy model
calculation to show the validity of this method. Finally, we show an interesting feature of this method
in view of the centrality resolution by using a toy model.

CHAPTER 4. VOLUME FLUCTUATION

53

Figure 4.2: κσ 2 (C4 /C2 ) of net-proton distribution as a function of centrality with different centrality
definition using UrQMD model [22].

4.2.1

Derivation

We assume that the independent particle productions of net-particle from the source of Npn . Thenr-th moment for each Npn can be expressed as

 r
d
r
h∆n i =M∆n (t)
,
(4.3)
dtr
t=0
where ∆n and M∆n denote the number of net-particle and the moment generating function for each
Npn . Due to the independent particle production model, a moment generating function for all source
will be expressed by the product of the moment generating functions from each source asM∆N (t) = [M∆n (t)]Npn ,(4.4)

where ∆N is the number of net-particle from all source ∆N = ∆n × Npn . Then r-th moment can beextracted by r derivatives on the moment generating function. The first and second order moments
are obtained asX
dM∆N (t)h∆N i =
P (Npn )
(4.5)
dt
t=0
Npn


X
Npn −1 dM∆n (t)=P (Npn ) Npn [M∆n (t)]= hNpn ih∆ni,
(4.6)
dt
t=0
Npn
 2

X
d M∆N (t)2h∆N i =
P (Npn )= hNpn (Npn − 1)ih∆ni2 + hNpn ih∆n2 i,
(4.7)
dt2
t=0
Npn

CHAPTER 4. VOLUME FLUCTUATION

54

where we apply additional summation over the probability distribution of participant nucleons P (Npn )to consider the fluctuating number of Npn . Once we obtain moments up to fourth order, they areconverted to cumulants:h∆N ic = hNpn ih∆ni,
2h∆N ic =h∆N 3 ic =4h∆N ic =

2

(4.8)
2

2
hNpn ih∆n ic + h∆ni hNpn
ic ,2
3
hNpn ih∆n3 ic + 3h∆nih∆n2 ic hNpn
ic + h∆ni3 hNpn
ic ,
4
3
2
2 2
2
hNpn ih∆n ic + 4h∆nih∆n ic hNpn ic + 3h∆n ic hNpn
ic34+6h∆ni2 h∆n2 ic hNpn
ic + h∆ni4 hNpn
ic .

(4.9)
(4.10)(4.11)

r i can be determined by someSince we can calculate h∆N r ic from the experimental data and hNpncmodels like Glauber, cumulants for each source h∆nr i are iteratively extracted from the first to r-thorder cumulant. Cumulants for all source are obtained ash∆N r ic = hNpn ih∆nr ic ,(4.12)

due to the assumption of independent source of participant nucleons. However, such derivation interms of cumulants for net-particles and participant nucleons becomes complicated and cumbersome
because of the conversions of explicit terms from moments to cumulants. In order to derive formulas
up to sixth order cumulant, thus we introduce factorial moments for participant nucleons fr andmoments for net-particles from each source µr = h∆nr i because terms related to participant nucleons
always appear as factorial moments as are seen in the first term in the right hand side in Eq. (4.7).
Measured moments up to sixth order can be expressed ash∆N i = f1 µ1 ,2h∆N i =
3h∆N i =h∆N 4 i =5h∆N i =
6h∆N i =

f2 µ21 + f1 µ2 ,f3 µ31 + 3f2 µ1 µ2 + f1 µ3 ,f4 µ41 + 6f3 µ21 µ2 + 3f2 µ22 + 4f2 µ1 µ3 + f1 µ4 ,f5 µ51 + 10f4 µ31 µ2 + 15f3 µ1 µ22 + 10f3 µ21 µ3 + 10f2 µ2 µ3 + 5f2 µ1 µ4 + f1 µ5 ,f6 µ61 + 15f5 µ41 µ2 + 45f4 µ21 µ22 + 15f3 µ32 + 20f4 µ31 µ3 + 60f3 µ1 µ2 µ3 + 10f2 µ23+15f3 µ21 µ4 + 15f2 µ2 µ4 + 6f2 µ1 µ5 + f1 µ6 .

(4.13)(4.14)(4.15)(4.16)(4.17)(4.18)

Moments from each source can be obtained up to sixth order by resolving Eqs.(4.13)–(4.18) iteratively.
After we convert them to cumulants from each source, corrected cumulants are obtained by Eq. (4.12).
Correction formulas more than fourth order cumulants are given byhih∆N 5 ic = hNpn iκ5 (∆n) + 5κ4 (∆n)κ1 (∆n) + 10κ3 (∆n)κ2 (∆n) κ2 (Npn )hi+ 10κ3 (∆n)κ21 (∆n) + 15κ22 (∆n)κ1 (∆n) κ3 (Npn ) + 10κ2 (∆n)κ31 (∆n)κ4 (Npn )h∆N 6 ic

+κ51 (∆n)κ5 (Npn ),(4.19)
hi= hNpn iκ6 (∆n) + 6κ5 (∆n)κ1 (∆n) + 15κ4 (∆n)κ2 (∆n) + 10κ23 (∆n) κ2 (Npn )hi+ 15κ4 (∆n)κ21 (∆n) + 60κ3 (∆n)κ2 (∆n)κ1 (∆n) + 15κ32 (∆n) κ3 (Npn )hi+ 20κ3 (∆n)κ31 (∆n) + 45κ22 (∆n)κ21 (∆n) κ4 (Npn ) + 15κ2 (∆n)κ41 (∆n)κ5 (Npn )+κ61 (∆n)κ6 (Npn ),which can be easily derived by using cumulant expansion technique [43].

(4.20)

CHAPTER 4. VOLUME FLUCTUATION

4.2.2

55

Toy model

The setup is following. We make the multiplicity distribution by using the Glauber and two component,and define the centrality using the multiplicity distribution. For proton productions, we tried two cases.
In one case, protons and antiprotons are generated from event by event Npart determined by Glaubermodel. In the other case, Npart is fixed to the averaged value in one centrality. Figure 4.3 shows the
net-proton distribution for two cases, and their ratios at each 0–5, 20–30 and 70–80% centralities. It
can be found that the shape of the net-proton distributions are affected by participant fluctuations.
Cumulants up to fourth order are shown in Fig. 4.4. For black circles Npart fluctuates which includes
participant fluctuation given by the Glauber model For red squares Npart is fixed, which are true valueof cumulant in this model. Black solid lines are results of VFC, and blue crosses are results of CBWC.
We see that VFC results gives the true value of cumulants up to fourth order. For CBWC, however,there are still some fractions of volume fluctuations compared to the true values for C3 and C4 .

Figure 4.3: Upper panel shows the net-proton distributions at 0–5, 20–30 and 60–70% centralities.
Npart fluctuates according to the Glauber model for red lines, while Npart is fixed in black lines. Lower
panels show the ratio between two net-proton distributions.

4.3

CBWC vs VFC

4.3.1

Centrality resolution effect

As discussed in Sec. 4.1, CBWC depends on the centrality resolution. In order to the see how VFC
behave with respect to the centrality resolution, we performed a simple toy model calculation. Although the setup for toy model is mostly the same with Sec. 4.2.2, three multiplicity distributions
are defined which have different centrality resolution as shown in Fig. 4.4, where we defined the red
multiplicity distribution as 100% centrality resolution, then multiplicity distributions which have 50%(blue) and 20% (green) centrality resolution were generated by random sampling from the multiplicity
distribution with 100% centrality resolution. We apply both CBWC and VFC in these three cases ofCHAPTER 4. VOLUME FLUCTUATION

56

Figure 4.4: Cumulants up to fourth order as a function of centrality. For black circles Npart fluctuates
which includes participant fluctuation given by the Glauber model For red squares Npart is fixed,
which are true value of cumulant in this model. Black solid lines are results of VFC, and blue crosses
are results of CBWC.
centrality resolution, then compare results between two methods and different centrality resolution.
Results are shown in the right hand side panels in Fig. 4.5. It can be found that the results of C3 andC4 without CBWC nor VFC and with CBWC are both enhanced by reducing the centrality resolution.
On the other hand, VFC results does not depend on the centrality resolution and give true value of
cumulants.

4.3.2

UrQMD model

Based on the discussions so far, let me summarize the merits and demerits for CBWC and VFC asfollowings. CBWC can suppress the volume fluctuation without any model assumption, but it depends
on the centrality resolution. On the other hand, VFC can eliminate the volume fluctuation, but it
requires some model assumptions. We start from the Glauber model and assume the independent
particle production model. In real experiments, however, the centrality resolution depends on thedetectors used for the centrality determination, and the independent particle production model could
be possibly broken, since lots of results from the heavy ion colliding experiments cannot be explainedby the superposition of p+p collisions. Therefore, both CBWC and VFC would never provide true
cumulants. Then the question is which correction will more closer to the true cumulants than theother. In order to obtain some hints for this question, the UrQMD data has been analyzed. The
UrQMD is an event generator that are tuned to reproduce single particle distributions measured
in heavy ion colliding experiments. The independent particle production model is expected to bebroken in UrQMD, therefore the comparison between CBWC and VFC will be better test compared
to previous models. In addition, we can define true cumulants by using the information on the eventby event Npart given by UrQMD. The centrality is defined by using pions and kaons in |η| < 1.0 as
shown in Fig. 4.6–(a). Protons and antiprotons in 0.4 < pT < 2.0 GeV/c and |y| < 0.5 are used forthe cumulant analysis. Figure 4.6–(b) shows the event by event Npart distribution given by UrQMD
for each centrality. The two-dimensional distributions in Np and Npbar at Npart = 240 are shown in
Fig. 4.6–(c) and (d) for 10-20 and 0-10% centrality, respectively. Then we can calculate cumulants
at each Npart bin in different centralities. Taking average of cumulants for each Npart bin in one

CHAPTER 4. VOLUME FLUCTUATION

57

Figure 4.5: (Left) Multiplicity distributions with 100% (red), 50% (blue) and 20% (green) centrality
resolution generated by Glauber and two component model. Distributions with 50% and 20% centrality
resolution were defined by random sampling from the distribution with 100% resolution. (Right)
Results up to fourth order cumulant as a function of centrality without CBWC nor CBWC, with
CBWC and with VFC from left to right. Dotted lines represent the true cumulants.

CHAPTER 4. VOLUME FLUCTUATION

58

centrality provides exactly the true cumulants in view of the participant fluctuation. This method
corresponds to the CBWC on Npart . In order to distinguish this method with the conventional CBWC
on the multiplicity distribution, I name both methods as CBWC-N and CBWC-M, respectively. The
UrQMD data has been analyzed with CBWC-N, CBWC-M and VFC. For VFC, the multiplicity
distribution in Fig. 4.6–(a) is fitted by the Glauber and two-component model to extract the Npart
distribution and its cumulants, as is done in the experiment, instead of using the event by event Npart
distribution given by UrQMD. The results from these three methods are shown in Fig. 4.7 for up tosixth order cumulant as a function of centrality. The cumulant ratios and the ratios of CBWC-M and
VFC to CBWC-N are shown in Fig. 4.8. It can be found that the both CBWC-M and VFC don’t
provide the true cumulants (CBWC-N) as expected. Therefore for the final analysis of C6 /C2 in theexperimental data, the both methods CBWC-M and VFC have been adopted.
In order to calculate the true cumulants defined as CBWC-N, we need to know event by event Npart .
One idea is to know the event by event Npart is to measure the spectator nucleons. But this is
impossible in the collider experiments since we can only measure a part of the fragments of spectator
nucleons by ZDC. On the other hand, however, it might be possible in the fixed-target experiment,
where all spectators can be measured by ZDC in principle. If the ZDC resolution is high enough tomeasure the energy of spectators with unit nucleons, the true cumulants will be achieved without any
assumptions.

Figure 4.6: (a) The refmult3 distribution and (b) the Npart distribution, where the colors represent thedifferent centralities. (c) and (d) show the two-dimensional histogram for Np and Npbar at Npart = 240
in 10-20 and 0-10% centralities.

CHAPTER 4. VOLUME FLUCTUATION

59

Figure 4.7: Cumulants up to sixth order as a function of centrality for CBWC-N (true), CBWC-M
and VFC.

Figure 4.8: C3 /C2 , C4 /C2 and C6 /C2 as a function of centrality for CBWC-N (true), CBWC-M and
VFC. The ratios of CBWC-M and VFC to CBWC-N are shown in bottom panels.

Chapter 5

Results and Discussions
5.1

Systematic uncertainties

The systematic studies has been done by changing the cut variables and the value of efficiency used
for the efficiency correction, and applying the unfolding with the beta-binomial model, which aresummarized in Tab. 5.1. First, the unfolding with beta-binomial model has been performed. As theunfolding method are incompatible with CBWC, the results of the efficiency correction and unfolding
are compared without any volume fluctuation corrections, which is shown in Fig 5.1. It is foundthat the unfolding results are consistent with the efficiency correction within statistical errors. The
non-binomial detector effects on C6 /C2 is thus negligible. The systematic uncertainties are estimatedby changing parameters for PID cuts, track quality cuts, pileup rejection, and the efficiency correction. Loosening the nσp and m2 cuts worsen purity of protons and antiprotons, and decreases kaon
contaminations. Requirement of small value of nHitsFit will increase the fraction of track splittings.
Large DCA cut value will increase secondary protons. Pileup events are included by loosening theslope in a plane of the number of tofmatched tracks and Refmult. The values used in the efficiency
correction are changed ±5% separately and simultaneously for low and high pT regions. Finally, thesystematic uncertainties are estimated by following formulas:
v

 #
"
u
sX
u 1 X C6 /C2 i,j − C6 /C2 def 2t
σsys = C6 /C2 def
Rj2 , Rj =
,
(5.1)
n
C
/C
6
2
def
j
i
C6
C6
where ( C)def represents the results with the default cuts, and ( C)i represents ith change of the cut
22
on jth variable.

60

CHAPTER 5. RESULTS AND DISCUSSIONS

VariablenσpDefault cut
< 2.0

m2

0.6 < m2 < 1.2

nHitsFit
DCA
nTofMatch =a × Refmult + b
nTofMatchBeta =c × Refmult + d

> 20
< 1.0
(a, b, c, d) =(0.5, −13, 0.46, −10)

Efficiency(εlowpT , εhighpT )Detector model

Binomial

61

Changed cut< 2.5
0.7 < m2 < 1.3
0.8 < m2 < 1.4
> 15
< 1.5
(0.25, −13, 0.46, −10)
(0.50, −13, 0.23, −10)(0, 0, 0, 0)(1.05 × εlowpT , 1.05 × εhighpT )(0.95 × εlowpT , 0.95 × εhighpT )(1.05 × εlowpT , 0.95 × εhighpT )(0.95 × εlowpT , 1.05 × εhighpT )Beta-binomial

effectworse purity
decreases
kaon contaminations
increases track splittings
increases secondary protons
increases pileup events

N/A
correct the non-binomial
detector efficiency

Table 5.1: Summary of the systematic checks.

Figure 5.1: Centrality dependence of C6 /C2 in Run11 without any volume fluctuation corrections.
Results from the efficiency correction are shown in red circles, and results from the unfolding with thebeta-binomial model are shown in green crosses.

CHAPTER 5. RESULTS AND DISCUSSIONS

5.2

62

Up to 4th order cumulants

Before the discussing the sixth order cumulant, let us compare results up to fourth order cumulant
between CBWC and VFC. Figure 5.2 shows the cumulants up to fourth order as a function of centrality
without CBWC nor VFC, with VFC and with CBWC. In the former two cases, results with different
centrality bin width of 8bin, 16bin and 32bin (corresponding to 2.5%, 5.0% and 10% centrality step)
are shown. It can be found that C3 and C4 in VFC results seems artificially suppressed at central
collisions with 8bin centrality class, while they become smooth and converge to certain value with small
centrality bin width of 16bin and 32bin. At first glance, it seems interesting because such centrality
dependence is not observed in a toy model, and one might expect some physics in experimental data.
However, it can be simply explained by multiplicity dependent efficiency. According to the analytical
calculations in Ref. [41], if the averaged efficiency is used under the situation of multiplicity dependent
efficiency, cumulants are artificially suppressed than they should be. One does know the fact that theefficiency depends on multiplicity as shown in Fig. 3.14. As was noted, the efficiency correction isapplied using the averaged efficiency in one centrality, then subtract the participant fluctuations. Thisis exactly the case discussed in Ref. [41], which indicates that the results of C3 and C4 without CBWC
or VFC shown in Fig. 5.2 are artificially suppressed. One can expect that such effect from multiplicity
dependent efficiency becomes small with narrow centrality bin width, because efficiency variations
will be reduced. Then the observed convergence of VFC results for C3 and C4 with 16bin and 32bin
centrality classes is consistent with the expectation, and it indicates that the effect of multiplicity
dependent efficiency becomes negligible with those centrality bin width.

Figure 5.2: Cumulants up to fourth order as a function of centrality using the experimental data.
Results without CBWC nor VFC, with VFC and with CBWC. In the former two cases, results withdifferent centrality bin width of 8bin, 16bin and 32bin (corresponding to 2.5%, 5.0% and 10% centralitystep) are shown.
Next, in order to see how much the STAR preliminary results on 3rd and 4th order fluctuations√
could be changed by using VFC instead of CBWC, datasets at sNN = 7.7, 19.6 and 62.4 GeV havebeen analyzed as well as 200 GeV. Figure 5.3 shows the centrality dependence of C3 /C2 normalized
by the value of Skellam and C4 /C2 at different beam energies. It can be found that VFC results of(C3 /C2 )/Skellam are systematically suppressed compared to CBWC in all centralities. The differences
between CBWC and VFC seem large in mid-central collisions, and small in peripheral and central

CHAPTER 5. RESULTS AND DISCUSSIONS

63

collisions, which is consistent with results of the toy model in Fig. 4.4. This trend of downward convex
in the centrality dependence changes to flat in high beam energy. For C4 /C2 , results are consistent√
between CBWC and VFC in central collisions at sNN = 7.7, 19.6 and 62.4 GeV. Next, let us discuss
the difference between the VFC results and the current STAR preliminary results. In Fig. 5.4, theVFC results are superimposed to the STAR preliminary results of beam energy dependence of C4 /C2
and Sσ/Skellam = (C3 /C2 )/Skellam in 0-5, 5-10 and 70-80% centralities [23]. It can be found that the
VFC results of κσ 2 are consistent with preliminary results within errors in 0-5 and 5-10% centralities.
The conclusion of non-monotonic behaviour in the beam energy dependence of κσ 2 would be thus
robust to the method for the volume fluctuation correction. However, large differences are observedin Sσ/Skellam and in 70-80% centrality of κσ 2 . Further understandings on the volume fluctuations
would be necessary to interpret these differences.

Figure 5.3: The centrality dependence of C3 /C2 normalized by the value of Skellam (top) and C4 /C2√(bottom) at sNN = 7.7, 19.6, 62.4 and 200 GeV. Results using CBWC are shown in blue crosses,
and VFC results are shown in green triangles. The Skellam baselines are shown in dotted lines.

5.3

6th order fluctuation

First, the Run10 and Run11 datasets are analyzed separately, and checked the consistency. Figure 5.5
shows the centrality dependence of C6 /C2 for Run10 and Run11, where the efficiency correction and
CBWC are applied. Due to the huge errors in central collisions, the expanded plot is shown in righthand side from 20% centrality. It can be found that they are mostly consistent within errors except
40-50% centrality that would be the statistical fluctuation. Then the both results are merged to reduce

CHAPTER 5. RESULTS AND DISCUSSIONS

64

Figure 5.4: The STAR preliminary results on the beam energy dependence of κσ 2 = C4 /C2 and√Sσ/Skellam = (C3 /C2 )/Skellam using CBWC [23], where VFC results at sNN = 7.7, 19.6 and 62.4
and 200 GeV are superimposed in green circles, orange diamonds and purple stars for 0-5, 5-10 and
70-80% centralities, respectively.

CHAPTER 5. RESULTS AND DISCUSSIONS

65

errors, which is shown in Fig. 5.6 with Skellam and the binomial baselines. It is found that the results
are systematically smaller than Skellam baseline, and mostly consistent with the binomial expectation.
It is natural for the binomial expectation can explain the results better than Skellam baseline since thebinomial distribution has two parameters while Skellam has one. One can see that the experimental
data is closer to the binomial distribution than Skellam distribution. In central collisions, however,results are still systematically suppressed compared to the binomial expectation. Since the results
contain the baryon stopping effect which is known to lead to the suppression with respect to the
Skellam baseline, one should compare the results with UrQMD model. CBWC is applied in UrQMD
model as is done in the experimental data, which is shown in Fig. 5.7 in purple band. It can befound that the results of UrQMD is systematically suppressed than the Skellam baseline, which canbe understood by the baryon stopping effect. But the experimental results still show smaller results
than UrQMD.

Figure 5.5: Centrality dependence of C6 /C2 for Run10 and Run11. Statistical errors are shown inbars, while systematic uncertainties are shown in brackets. The expanded plot is shown in right handside from 20% centrality.

Figure 5.6: C6 /C2 as a function of centrality, where results from Run10 and Run11 are merged. The
expanded plot is shown in right hand side from 20% centrality. Skellam and the binomial baselines
are also shown.

CHAPTER 5. RESULTS AND DISCUSSIONS

66

Figure 5.7: C6 /C2 as a function of centrality with Skellam baseline and the results of UrQMD model.
CBWC is applied in UrQMD model as is done in experimental data.

Figure 5.8: C6 /C2 as a function of centrality using CBWC or VFC for both experimental data and
UrQMD model. Results with CBWC are shown in purple and the results with VFC are in green.
Results of UrQMD model are shown in band.

CHAPTER 5. RESULTS AND DISCUSSIONS

5.4

67

Interpretation

When one discusses the centrality dependence of higher order fluctuations, the finite size effect should
be considered. In heavy ion collisions the system size is limited in the size of nucleus, which indicatesthat the fluctuations don’t substantially grow up. In other words, number of produced particles issmall in peripheral collisions, thus anything other than statistical fluctuations would not be observed,
while the fluctuations will get larger in central collisions (but still restricted to the size of nuclei).
Another issue is that one still doesn’t know from which centrality QGP is created. If there is no QGP,
off course, phase transition should not occur, which leads to null fluctuation. In order to see whether
some signals from the phase transition are observed or not, the significance with respect to zero isplotted as a function of centrality in Fig. 5.9. For CBWC results, it is found that the C6 /C2 shows thesystematic reduction from positive to negative values in peripheral and semi-central collisions. It canbe understood by the absence of QGP or finite size effect. Even in semi-central collisions, QGP might
be created and C6 /C2 would become more sensitive to the phase transition than finite size effect, then
the negative values of C6 /C2 are observed.

Figure 5.9: Significance with respect to zero as a function of centrality for CBWC and VFC. The
expanded plot is shown in right hand side panel.

5.5

Comparison with theoretical calculations

The model calculation in Ref. [14] would provide us hints to restrict the location of the chiral crossoverwith respect to the chemical freeze-out. The results of C4 /C2 and C6 /C2 on net-proton distribution
B
B
B
in 0-10 and 30-40 % centralities, as well as a rough estimates on χB
6 /χ2 and χ4 /χ2 with different
freeze-out conditions from model calculation [14] are summarized in Tabs. 5.2 and 5.3. It can be foundthat the experimental results are consistent with the freeze-out condition of T freeze /Tpc ' 1. This isalso consistent with the conjecture in Sec. 1.2.3 that the crossover region might be close enough tothe chemical freeze-out line.
One should note that the absolute values in Fig. 1.9 is model dependent, quantitative comparisonsB
can thus only be done with the lattice QCD calculation. Recently the value of χB
6 /χ2 has been
calculated by the lattice QCD [24], which is shown in Fig. 5.10 as a function of the temperature. At√
sNN , the chemical freeze-out temperature is around 166 MeV. Corresponding results from the lattice
QCD are superimposed to the experimental results as red shadow in Fig. 5.11. It can be found thatthe experimental results do not conflict with the lattice QCD, but there are still large errors for bothresults. It should be also noted that the direct comparison with experiments and the lattice QCD

CHAPTER 5. RESULTS AND DISCUSSIONS
Centrality
C4 /C2
C6 /C20-10 %
0.98 ± 0.11
0.76 ± 0.11
−0.97 ± 18.4
2.98 ± 12.0

30-40 %
0.74 ± 0.03
0.42 ± 0.01
−2.64 ± 1.63
−3.57 ± 0.95

Table 5.2: Results on C4 /C2 and C6 /C2 in
0-10 and 30-40 % centralities at Au+Au col√lisions sNN = 200 GeV. CBWC and VFC
results are shown in upper and lower at eachcell. (Statistical errors only)

68

B
χB
4 /χ2
B
χ6 /χB2

T freeze /Tpc ' 1
∼ 0.5
<0

T freeze /Tpc ≤ 0.9≥1≥1

Table 5.3: Rough estimates on the ratio offourth to second and sixth to second orderbaryon number susceptibilities with different
freeze-out conditions.

would be difficult. There are several experimental effects that are not fully considered in the theory,
such as the global baryon number conservation effect [42, 44], volume fluctuation effects [22, 42] andthe difference between net-baryon with net-proton [36].

B
Figure 5.10: χB
the lattice QCD as a function of temperature [24]. The red shaded
6 /χ2 calculated with√area it the region of interest at sNN = 200 GeV.

CHAPTER 5. RESULTS AND DISCUSSIONS

69√
Figure 5.11: Final results on C6 /C2 of net-proton distribution at sNN = 200 GeV, where thecorresponding results calculated by the lattice QCD on the ratio of sixth to second order susceptibilityof the net-baryon distribution [24] are superimposed in red shadows.

Chapter 6

Conclusions
The ratio of the sixth to second order cumulant of net-proton multiplicity distributions in Au+Au√
collisions at sNN = 200 GeV at the STAR experiment has been measured as a function of centrality. The experimental non-binomial detector efficiency was determined by using the embedding
simulations. The non-binomial model was implemented into the unfolding, and it was found thatthe non-binomial detector effect was negligible, which indicates that the analytical correction for thebinomial detector efficiency is applicable to the data. Regarding the volume fluctuation effect arisingfrom the participant fluctuations, it was found that both CBWC and VFC would not be correct by√using the UrQMD model. VFC was applied to C3 /C2 and C4 /C2 at sNN = 7.7, 19.6, 62,4 and
200 GeV. In most centralities and beam energies, VFC results are systematically smaller than CBWC
results. It was found that the non-monotonic behaviour, which could be possible signature of the critical end point, observed in C4 /C2 as a function of beam energy using CBWC was robust even if VFC
is applied instead of CBWC. Both CBWC and VFC were applied to C6 /C2 . For CBWC, the results
show positive values in the peripheral collisions and negative values in central collisions systematically.
This could be possibly a signature of the crossover phase transition in central collisions. For VFC,on the other hand, the results show negative values for all centralities, which could indicate that the
QGP is formed and the phase transition is probed even in peripheral collisions. The results show
systematically smaller values compared to the statistical baselines and the UrQMD model. Based on
the theoretical calculation, observed negative values of C6 /C2 would indicate that the chiral crossover
region is located near the chemical freeze-out line. The result was also compared with the lattice QCD.
The experimental and lattice QCD results are consistent with large errors. However, one should noteseveral effects like volume fluctuations, global baryon number conservation, and the difference between
the net-baryon and net-proton, which are not fully understood. Further development of experiments
and theories would be essential to extract more definite information on the QCD phase structure.

70

Appendix A

Centrality definition
Based on Sec. 4.1, we used charged pions and kaons measured by TPC in |η| < 1.0, that we call
refmult3, for the centrality determination. Figure A.1 shows averaged refmult3 as a function of
ZDC conincidence rate which represents the luminosity. We can see that refmult3 decreases withincreasing the luminosity as shown in blue markers, which was corrected as shown in red markers.
After luminosity correction, we looked at the refmut3 distributions as shown in the left hand side in
Fig. A.2 to perform fitting with following function:f (x) = a × Erf(−b × (x − c)) + a,Z x
22
Erf(x) = √e−t dt,
π 0

(A.1)
(A.2)

where Erf(x) is the error function and c was extracted as a parameter which represents the maximum
value of refmult3 distribution. This was performed at different Vz windows with 2 cm step, which isplotted as a function of Vz in the right hand side in Fig. A.2, where we observe Vz dependence which
was corrected by 2nd polynomial fitting.

Figure A.1: Luminosity dependence of averaged refmult3 in −2 < V z < 0 cm. Blue markers areuncorrected and fitted by first order polynomial function (red lines) which is used for the correction.
Corrected results are shown in blue markers.
After luminosity and acceptance corrections, we define the centrality. Since trigger efficiency drops
in peripheral collisions, experimental measured multiplicity distribution cannot be directly used to di71

APPENDIX A. CENTRALITY DEFINITION

72

Figure A.2: (Left) Tail of the refmult3 distribution fitted by error function shown in a red line.
The blue dotted line represents the parameter c in Eq. (A.1). (Right) Maximum value of refmult3
distribution as a function of Vz (black) and the corrected distribution (red).
vide the centrality. So we have to perform some model calculations to reproduce the experimental
measured multiplicity distribution except very peripheral collisions, which is used to divide the centrality. We used Glauber and two-component model to generate the multiplicity distributions. In Glauber√model, inelastic cross section for p+p collisions are set to σpp = 43.6 (mb) at sN N = 200 GeV.
Other parameters for geometry are set to be of Au nuclei. 100 M events were proceeded. Once thenumber of participant nucleons (Npart ) and the number of binary collisions (Ncoll ) are determined by
Glauber model, particles were generated by two-component model with negative binomial fluctuation:


Npart
Nprod = npp (1 − x)
+ xNcoll ,(A.3)2(µ/k)aΓ(a + k)
(A.4)
npp =
Γ(a + 1)Γ(0) (1 + µ/k)−a−kwhere x is the fraction of hard component in particle production, npp represents the number of
produced particles per unit p+p collision, and Γ represents the gamma function. In each collision
event, the number of ”source” for particle production is determined according to the second terms
shown in the first equation, then particles are generated from the ”source” by npp , which is randomly
selected by the negative binomial distribution in the second equation. Multiplicity dependent efficiency
were also implemented to reproduce the multiplicity distribution at the STAR experiment:
Nmeas = εmult × Nprod ,


Nprod
εmult (Nprod ) = 0.98 1 − εconst
,
560

(A.5)
(A.6)

where the slope of Nmeas distribution can be controlled by εconst . Those parameters were selected sothat the Nmeas distributions can well reproduced the refmult3 distribution in refmult3 > 100, whichis shown in the left hand side in Fig. A.3, where red points represent Nmeas distribution and black
points are refmult3 distribution. Right hand side panel in Fig. A.3 shows the ratio of Nmeas to refmult3
distributions. We can see the fitting is good (χ2 /ndf ) except refmult3 < 100. The blue line is thefitting in refmult3 < 100, which can be used for trigger inefficiency correction for cumulant analysis.
MC multiplicity distribution are divided into 9 centrality classes, 0-5%, 5-10% and 10% step up toAPPENDIX A. CENTRALITY DEFINITION

73

70-80% centrality. Table A.1 shows the results of centrality definition, where low edge multiplicity bin
and the averaged number of participant nucleons for each centrality class are shown.

Figure A.3: (Left) Refmult3 distribution (black points) fitted by MC distribution generated by Glauber
and two-component model (red points). (Right) Ratio of MC to refmult3 distribution. The solid blueline represents the fitting for trigger inefficiency correction.

Centrality(%)
0–5
5–10
10–20
20–30
30–40
40–50
50–60
60–70
70–80

Low edge bin
784
674
484
332
216
131
73
36
15

hNpart i
350
300
235
167
115
75
47
26
13

√
Table A.1: Summary for the centrality definition at sN N = 200 GeV in run11 datasets. The middle
row shows the low edge refmult3 bin and the right row is the averaged number of participant nucleons.

Appendix B

Tips on the efficiency correction
B.1

Stirling number of the first kind

The falling factorial moment is defined by using Stirling number of the first kind as:(x)n =nXs(n, k)xk ,

(B.1)

k=0

with s(n, k) being the Stirling number of the first kind, which is given by   nn−k n
= |s(n, k)|.
,s(n, k) = (−1)
k
k

(B.2)

The second equation represent the unsigned Stirling number of the first kind, which is calculated by
the following recursive relation:


  

n+1
n
n
=n
+
,
(B.3)
k
k
k−1
with the initial condition

B.2

 
   
0
n
0
= 1,
=
= 0.
0
0
n

(B.4)

Factorial moments with the binomial model

In this section, I introduce Eq. (3.5) which holds under the binomial model. For convenience, I use
the Pochhammer symbol(x)i = x(x − 1)(x − 2) · · · (x − i + 1) =

74

x!
.
(x − i)!

(B.5)

APPENDIX B. TIPS ON THE EFFICIENCY CORRECTION

75

Starting from Eq. (3.3),
fab =

X
a

s(a, i)M

i

i=0

=

X

b
X

s(b, j)M̄

j



j=0

p(M, M̄ ) × (M )a (M̄ )b

M,M̄

=

X

(M )a (M̄ )b

M,M̄

=

X

P (N, N̄ )

N,N̄

=

=

X

X

P (N, N̄ )

N,N̄

X
M,M̄

M!
N!
M̄ !
N̄ !
εM ε̄M̄ (1 − ε)N −M (1 − ε̄)N̄ −M̄
(M − a)! M̄ − b M !(N − M )! M̄ !(N̄ − M̄ )P (N, N̄ )(N )a (N̄ )b

X

N,N̄

M,M̄

X

a b

P (N, N̄ )(N )a (N̄ )b ε ε̄N,N̄

= εa ε̄b Fab .

N̄ !
N!
ε̄M̄ (1 − ε̄)N̄ −M̄
εM (1 − ε)N −M
M !(N − M )!
M̄ !(N̄ − M̄ )!

(N − a)!
(N̄ − b)!
εM ε̄M̄ (1 − ε)N −M (1 − ε̄)N̄ −M̄
(M − a)!(N − M )! (M̄ − b)!(N̄ − M̄ )!
∞
X
M =a

B(N − a, M − a; ε)∞X

B(N̄ − b, M̄ − b; ε̄)

M̄ =b

(B.6)

Appendix C

Net-proton distribution
In this appendix, the shape of the net-proton distribution is discussed. Corrected cumulants deviate
from the Skellam expectation. Furthermore, higher order cumulants more than second order have been
observed, which indicates that the observed net-proton distribution is no longer the Gaussian. The
unfolded distribution is compared with the Skellam and Gaussian distribution and discuss how thevalues of higher order cumulants up to the sixth order could appear as the shape of the distribution.

C.1

Shape for the sixth order cumulant

As was explained in Sec. 1.3, the third order cumulant represents the asymmetry of the distribution,and the fourth order cumulant represents the sharpness of the distribution. What about for the fifthand sixth order cumulant? For our intuitive understanding, let us consider to stack n boxes for then-th order cumulant, and suppose that the structure made by n boxes represents the shape of the
distribution. Figure C.1 shows a simple sketch for up to the sixth order. For one box, the box can
only be located there without any other choices. For two boxes, two boxed can be placed vertically or
horizontally, the entire shape of two boxes can thus represent the width of the distribution. By adding
another box, for three boxes, the last one box can be placed at the right or left hand side of verticallystacked two boxes, which represents the asymmetry of the distribution. For four boxes, box-like andconvex upward shape can be considered. In the case of five and six boxes, the degree of freedom at
the top and side part of the distribution is added with respect to the cases of three and four boxes,
respectively. Therefore, the ”expanded” or ”dented” shape at the side and top part of the distribution
could be possibly the reason for the finite value of the sixth order cumulant.
A simple toy model has been performed to check this conjecture. First, a Gaussian distribution
with mean and sigma being 0 and 5 is generated. In order to change the shape around the side
part of the Gaussian, additional two Gaussians are also prepared, which have the mean value of±15 and the sigma of 1. The original Gaussian is then modified to be expanded or dented shapeby adding or subtracting two symmetric Gaussians, which are shown in Fig. C.2. The ratio of themodified distribution to the Gaussian fitting in bottom panels. The kink structure is observed aroundN = 15, which represents the dented and expanded shape. Cumulants up to the sixth order on the two
distributions are shown in Tab. C.1. For the ordinary Gaussian distribution, higher order cumulants
become zero by definition. Even for the expanded and dented distributions, they are still symmetric,C3 and C5 are thus zero. But the even order cumulants are affected. For the expanded distribution,
the shape around the top becomes relatively sharper than the side part, so the value of C4 are positive.
On the other hand, the value of C4 of the dented distribution is negative. It can be also found thatthe value of C6 is negative for the expanded distribution, while positive for the dented one. Therefore,
76

APPENDIX C. NET-PROTON DISTRIBUTION

77

Figure C.1: A simple sketch to stack n boxes for the n-th order cumulant.
the observed negative value of C6 could be observed as the expanded shape around the side part ofthe distribution.
Cumulants
C1
C2
C3
C4
C5
C6

Gaussian
0
25
0
0
0
0

expanded
025.88 ± 0.01
077.39 ± 1.03
0−5890 ± 133

dented024.69 ± 0.01
0−40.17 ± 0.99
02131 ± 127

Table C.1: Cumulants up to the sixth order on the distributions in Fig. C.2.

C.2

Net-proton distribution

Let us compare the observed net-proton distribution with the conjectures discussed in previous section.
Top panels in Fig. C.3 shows the net-proton distributions in 0-5, 30-40 and 70-80 % centralities. Bluesquares are the measured distribution. Red circles are unfolded distribution assuming the binomial
model. Black lines represent Skellam distributions, and dotted green lines are Gaussian distributions.
The ratios of the unfolded distribution to the Skellam and Gaussian distributions are shown in thebottom panels. Since the unfolded distribution includes the effect of the volume fluctuation, corresponding cumulants are red open squares in Fig. 5.2 in which any volume fluctuation corrections are
not applied. Unfolded distributions are a bit narrower than the Skellam distributions in all centralities,
which leads to the shape of convex upward as is seen in the ratio in Fig. C.3. This is consistent withthe results on C2 . In 0-5 and 30-40% centralities, the ratio of the unfolded to the Skellam distribution
are skew to the left hand side. This asymmetry is also observed in the ratio of the unfolded to theGaussian distribution. Larger C3 than the Skellam baseline in Fig. 5.2 is thus observed. In addition,
the ratio of the unfolded to the Gaussian distribution behaves like ”W”. This indicates that theAPPENDIX C. NET-PROTON DISTRIBUTION

78

Figure C.2: Gaussian distributions with the mean value of 0 and the sigma of 5, in which another twosymmetric Gaussians around N = ±15 are added (left) or subtracted (right). Red lines represent theGaussian fittings. The ratio to the fitting are shown in the bottom panels.

APPENDIX C. NET-PROTON DISTRIBUTION

79

unfolded distribution is sharper than the Gaussian, which can explain the positive C4 . Interestingly,
some modifications like convex upward are observed at N = −5 and 20 in 0-5 %, and at N = −8 and
10 in 30-40 % centrality in the ratio of the unfolded to the Skellam distributions. Such modifications
cannot be explained by cumulants up to the fourth order, which could be the origin of the finite C6 in
30-40% centrality. One should note that the error bar of the unfolded distribution is not correct. The
current unfolding technique cannot allow us to calculate errors on the distribution but errors on the
cumulants. Thus the observed modification of the unfolded distribution might be negligible includingtrue errors.

Figure C.3: (Top) Net-proton distributions in 0-5, 30-40 and 70-80% centralities. Blue squares are themeasured distribution. Red circles are unfolded distribution without any volume fluctuation corrections. Black lines represent Skellam distributions, and dotted green lines are Gaussian distributions.
(Bottom) Ratios of the unfolded distribution to the Skellam and Gaussian distributions.

Appendix D

Robustness to the time dependence
According to Fig. 2.12, it can be found that the total multiplicity has time dependence, which is around±5 % in 3 σ. Assuming that the single particle efficiency directly depends on the total multiplicity, theefficiency variation could be assumed to be within ±5 %. This efficiency variation has been included in
the systematic studies. Here some hints will be given on how much robust the higher order cumulants
are, assuming that the efficiency drops in some fractions of total statistics by using a simple toy model.
Protons and antiprotons are generated with 100 M events with mean values of 5.64 and 4.39. The√
efficiencies are assumed to be ε = 0.65, which corresponds to 30-40 % centrality at sNN 200 GeV.
The efficiency is dropped with ∆ε in some ”bad” events, then apply the efficiency correction by using
ε = 0.65. Figure D.1 shows the relative difference to the true value of cumulants up to the sixth order
as a function of the fraction of bad events which have the efficiency of ε(1 − ∆ε) with ∆ε = 0.05,
0.2, and 0.4. It can be found that the effect is smaller than 1 % for ∆ε = 0.05 for all the order ofcumulants. For the case of ∆ε = 20 %, C3 and C4 start to deviate significantly from 5% fractionof bad events, whereas the effects on C5 and C6 are statistically negligible. As a simple criterion,
therefore, the fraction of bad events in which the efficiency drops 20 % has to be controlled to be lessthan 5 % in order to measure cumulants up to the fourth order with 1 % accuracy.

Figure D.1: Relative deviations with respect to the true value of cumulants up to the sixth order asa function of the fraction of bad events in which the efficiency drops with ∆ε.

80

Bibliography[1] Particle Data Group Collaboration, J. Beringer et al., Review of Particle Physics (RPP),
Phys. Rev. D86 (2012) 010001.
[2] CP-PACS Collaboration, A. Ali Khan et al., Light hadron spectroscopy with two flavors ofdynamical quarks on the lattice, Phys. Rev. D65 (2002) 054505, [hep-lat/0105015]. [Erratum:
Phys. Rev.D67,059901(2003)].
[3] K. Fukushima and T. Hatsuda, The phase diagram of dense QCD, Rept. Prog. Phys. 74 (2011)014001, [1005.4814].
[4] C. Nonaka and M. Asakawa, Modeling a Realistic Dynamical Model for High Energy Heavy Ion
Collisions, PTEP 2012 (2012) 01A208, [1204.4795].
[5] STAR Collaboration, L. Adamczyk et al., Bulk Properties of the Medium Produced in
Relativistic Heavy-Ion Collisions from the Beam Energy Scan Program, Phys. Rev. C96 (2017),
no. 4 044904, [1701.07065].
[6] STAR Collaboration, L. Adamczyk et al., Beam Energy Dependence of Jet-Quenching Effects√
in Au+Au Collisions at sNN = 7.7, 11.5, 14.5, 19.6, 27, 39, and 62.4 GeV, 1707.01988.
[7] STAR Collaboration, L. Adamczyk et al., Beam-Energy Dependence of the Directed Flow ofProtons, Antiprotons, and Pions in Au+Au Collisions, Phys. Rev. Lett. 112 (2014), no. 16
162301, [1401.3043].
[8] STAR Collaboration, L. Adamczyk et al., Beam energy dependence of moments of the
net-charge multiplicity distributions in Au+Au collisions at RHIC, Phys. Rev. Lett. 113 (2014)
092301, [1402.1558].
[9] STAR Collaboration, L. Adamczyk et al., Energy Dependence of Moments of Net-proton
Multiplicity Distributions at RHIC, Phys. Rev. Lett. 112 (2014) 032302, [1309.5681].
[10] STAR Collaboration, L. Adamczyk et al., Collision Energy Dependence of Moments of
Net-Kaon Multiplicity Distributions at RHIC, 1709.00773.
[11] X. Luo, Exploring the QCD Phase Structure with Beam Energy Scan in Heavy-ion Collisions,
Nucl. Phys. A956 (2016) 75–82, [1512.09215].
[12] M. A. Stephanov, On the sign of kurtosis near the QCD critical point, Phys. Rev. Lett. 107(2011) 052301, [1104.1627].
[13] S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo,
The QCD equation of state with dynamical quarks, JHEP 11 (2010) 077, [1007.2580].
81

BIBLIOGRAPHY

82[14] B. Friman, F. Karsch, K. Redlich, and V. Skokov, Fluctuations as probe of the QCD phasetransition and freeze-out in heavy ion collisions at LHC and RHIC, Eur. Phys. J. C71 (2011)
1694, [1103.3511].
[15] M. Asakawa and M. Kitazawa, Fluctuations of conserved charges in relativistic heavy ion
collisions: An introduction, Prog. Part. Nucl. Phys. 90 (2016) 299–342, [1512.05038].
[16] M. Anderson et al., The Star time projection chamber: A Unique tool for studying highmultiplicity events at RHIC, Nucl. Instrum. Meth. A499 (2003) 659–678, [nucl-ex/0301015].
[17] X. Zhang, Lifetime Calculations for Unstable Nuclei and φ Meson Production at RHIC. PhD
thesis, Nanjing University, 2010.
[18] W. J. Llope, The large-area time-of-flight upgrade for STAR, Nucl. Instrum. Meth. B241 (2005)
306–310.
[19] W. J. Llope et al., The STAR Vertex Position Detector, Nucl. Instrum. Meth. A759 (2014)23–28, [1403.6855].
[20] X. Luo, Search for experimental signatures of critical point with beam energy scan at STAR, in
CiRfSE workshop 2016.
[21] T. Nonaka, M. Kitazawa, and S. Esumi, More efficient formulas for efficiency correction of
cumulants and effect of using averaged efficiency, Phys. Rev. C95 (2017), no. 6 064912,[1702.07106].
[22] X. Luo, J. Xu, B. Mohanty, and N. Xu, Volume fluctuation and auto-correlation effects in the
moment analysis of net-proton multiplicity distributions in heavy-ion collisions, J. Phys. G40(2013) 105104, [1302.2332].
[23] X. Luo, Energy Dependence of Moments of Net-Proton and Net-Charge Multiplicity
Distributions at STAR, vol. CPOD2014, 2015. 1503.02558.
[24] A. Bazavov et al., The QCD Equation of State to O(µ6B ) from Lattice QCD, Phys. Rev. D95(2017), no. 5 054504, [1701.04325].
[25] M. Cheng et al., The Transition temperature in QCD, Phys. Rev. D74 (2006) 054507,[hep-lat/0608013].
[26] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, The Order of the quantum
chromodynamics transition predicted by the standard model of particle physics, Nature 443(2006) 675–678, [hep-lat/0611014].
[27] M. A. Stephanov, QCD phase diagram: An Overview, PoS LAT2006 (2006) 024,
[hep-lat/0701002].
[28] O. Philipsen, Status of the QCD Phase Diagram from Lattice Calculations, Acta Phys. Polon.
Supp. 5 (2012) 825–835, [1111.5370].
[29] STAR Collaboration, L. Adamczyk et al., Beam-energy dependence of charge separation along
the magnetic field in Au+Au collisions at RHIC, Phys. Rev. Lett. 113 (2014) 052302,[1404.1433].

BIBLIOGRAPHY

83[30] STAR Collaboration, L. Adamczyk et al., Beam Energy Dependence of the Third Harmonic of
Azimuthal Correlations in Au+Au Collisions at RHIC, Phys. Rev. Lett. 116 (2016), no. 11112302, [1601.01999].
[31] H. Stoecker, Collective flow signals the quark gluon plasma, Nucl. Phys. A750 (2005) 121–147,[nucl-th/0406018].
[32] R. V. Gavai and S. Gupta, Lattice QCD predictions for shapes of event distributions along thefreezeout curve in heavy-ion collisions, Phys. Lett. B696 (2011) 459–463, [1001.3796].
[33] A. Bzdak and V. Koch, Acceptance corrections to net baryon and net charge cumulants, Phys.
Rev. C86 (2012) 044904, [1206.4286].
[34] A. Bzdak and V. Koch, Local Efficiency Corrections to Higher Order Cumulants, Phys. Rev.
C91 (2015), no. 2 027901, [1312.4574].
[35] X. Luo, Unified Description of Efficiency Correction and Error Estimation for Moments ofConserved Quantities in Heavy-Ion Collisions, Phys. Rev. C91 (2015), no. 3 034907,[1410.3914].
[36] M. Kitazawa and M. Asakawa, Relation between baryon number fluctuations and experimentallyobserved proton number fluctuations in relativistic heavy ion collisions, Phys. Rev. C86 (2012)
024904, [1205.3292]. [Erratum: Phys. Rev.C86,069902(2012)].
[37] X. Luo, Error Estimation for Moments Analysis in Heavy Ion Collision Experiment, J. Phys.
G39 (2012) 025008, [1109.0593].
[38] A. Bzdak, R. Holzmann, and V. Koch, Multiplicity dependent and non-binomial efficiency
corrections for particle number cumulantsMultiplicity-dependent and nonbinomial efficiency
corrections for particle number cumulants, Phys. Rev. C94 (2016), no. 6 064907, [1603.09057].
[39] T. Adye, Unfolding algorithms and tests using RooUnfold, in Proceedings, PHYSTAT 2011
Workshop on Statistical Issues Related to Discovery Claims in Search Experiments andUnfolding, CERN,Geneva, Switzerland 17-20 January 2011, (Geneva), pp. 313–318, CERN,
CERN, 2011. 1105.1160.
[40] STAR Collaboration, B. I. Abelev et al., Identified baryon and meson distributions at large√transverse momenta from Au+Au collisions at sNN = 200 GeV, Phys. Rev. Lett. 97 (2006)
152301, [nucl-ex/0606003].
[41] A. Bzdak, R. Holzmann, and V. Koch, Multiplicity dependent and non-binomial efficiency
corrections for particle number cumulantsMultiplicity-dependent and nonbinomial efficiency
corrections for particle number cumulants, Phys. Rev. C94 (2016), no. 6 064907, [1603.09057].
[42] P. Braun-Munzinger, A. Rustamov, and J. Stachel, Bridging the gap between event-by-event
fluctuation measurements and theory predictions in relativistic nuclear collisions, Nucl. Phys.
A960 (2017) 114–130, [1612.00702].
[43] M. Kitazawa, Efficient formulas for efficiency correction of cumulants, Phys. Rev. C93 (2016),
no. 4 044911, [1602.01234].
[44] A. Bzdak, V. Koch, and V. Skokov, Baryon number conservation and the cumulants of the net
proton distribution, Phys. Rev. C87 (2013), no. 1 014901, [1203.4529].