A general procedure for detector-response correction of higher order cumulants
JJ-PARC-TH-0118
A general procedure for detector-response correction of higher order cumulants
Toshihiro Nonaka,
1, 2, ∗ Masakiyo Kitazawa,
3, 4, † and ShinIchi Esumi ‡ Key Laboratory of Quark & Lepton Physics (MOE) and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China Tomonaga Center for the History of the Universe, University of Tsukuba, Tsukuba, Ibaraki 305, Japan Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan J-PARC Branch, KEK Theory Center, Institute of Particle and Nuclear Studies,KEK, 203-1, Shirakata, Tokai, Ibaraki, 319-1106, Japan (Dated: May 9, 2018)We propose a general procedure for the detector-response correction (including efficiency correc-tion) of higher order cumulants observed by the event-by-event analysis in heavy-ion collisions. Thismethod makes use of the moments of the response matrix characterizing the property of a detector,and is applicable to a wide variety of response matrices such as those having non-binomial responsesand including the effects of ghost tracks. A procedure to carry out the detector-response correctionof realistic detectors is discussed. In test analyses, we show that this method can successfully re-construct the cumulants of true distribution for various response matrices including the one havingmultiplicity-dependent efficiency.
I. INTRODUCTION
Fluctuations are important observables in relativisticheavy-ion collisions for the search for the QCD criti-cal point and the phase transition to the deconfinedmedium [1, 2]. In particular, non-Gaussianity of fluc-tuations characterized by the higher order cumulants isbelieved to be sensitive to these phenomena [3–5]. Activemeasurements of fluctuation observables have been per-formed by the event-by-event analyses at RHIC [6–9], theLHC [10, 11], and the NA61 [12] and HADES [13] Col-laborations. Future experimental facilities, J-PARC [14],FAIR [15], and NICA [16], will also contribute this sub-ject.In relativistic heavy-ion collisions, the cumulants of aparticle-number distribution are obtained from an event-by-event histogram of the particle number observed ex-perimentally. However, because of the imperfect capa-bility of detectors the experimentally-observed event-by-event histogram is modified from the true distribution,and accordingly their cumulants are also altered by theseartificial effects. In the experimental analysis, this effecthas to be corrected. In this paper we call this procedureas the detector-response correction.Compared to standard observables given by expecta-tion values, the detector-response correction of the cu-mulants higher than the first order is more involved, be-cause the change of a distribution function modifies itscumulants in a non-trivial way [1]. So far, the detector-response correction has been discussed by focusing on theefficiency correction, i.e. the correction for the effects ofthe loss of particles at the measurement. It has been es-tablished that the correction can be carried out if oneassumes that the probability that the detector observesa particle (efficiency) is uncorrelated for individual parti-cles. In this case, the detector’s response is described bythe binomial distribution [1, 17], that we call the bino-mial model [17–24]. Recently, using the method proposed ∗ [email protected] † [email protected] ‡ [email protected] in Ref. [24], the analysis of the net-proton number cumu-lants is realized up to sixth order [25].The assumption for the binomial model, however, ismore or less violated at typical detectors in heavy-ioncollisions. First, the multiplicity dependence of the effi-ciency of realistic detectors [9, 13, 26] suggests the exis-tence of the correlations between efficiencies of differentparticles [27]. Even worse, typical detectors sometimesmeasure ghost tracks, i.e. non-existing particles. Theestimate on systematic uncertainty arising from these ef-fects is important for reliable experimental analyses ofhigher order cumulants.One of the general procedures for the detector-responsecorrection is the unfolding method [28–30]. This methodanalyzes the true distribution function from experimentalresults and knowledge on the detector. Strictly speaking,however, the reconstruction of the true distribution func-tion is an ill-posed problem. The estimate on the system-atic uncertainty of the final results is a nontrivial task inthis method, and a large numerical cost is required forthe iterative analysis to obtain the distribution. Further-more, the analysis of the distribution itself seems redun-dant because the cumulants are relevant quantities formany purposes. It is desirable to have a method whichenables the analysis of the cumulants directly withoutusing the distribution.In the present study, we propose a new method to per-form the detector-response correction of cumulants di-rectly. In this method, we relate the moments of true andobserved distributions without using the distribution ex-plicitly. This method can solve the detector-response cor-rection exactly for a wide variety of detector’s responses;for example, those parametrized by the hypergeomet-ric distribution and the binomial model with fluctuat-ing probability. By introducing an approximation with atruncation this method can also deal with the correctionof realistic detectors whose response is estimated only byMonte Carlo simulations. We demonstrate by explicit nu-merical analyses that the correction by this method canbe carried out successfully for non-binomial response ma-trices and the response matrix representing multiplicity-dependent efficiency.This paper is organized as follows. In Sec. II, we dis-cuss a general procedure for the detector-response correc- a r X i v : . [ phy s i c s . d a t a - a n ] M a y tion. The application of this method to the correction ofrealistic detectors is then discussed in Sec. III. We per-form test analyses of this method in Secs. IV and V: Wedeal with the exactly solvable cases in Sec. IV, and themultiplicity-dependent efficiency in Sec. V, respectively.The last section is devoted to a summary. II. DETECTOR-RESPONSE CORRECTIONA. Problem
Let us first clarify the problem. We consider a measure-ment of the cumulants of an event-by-event distributionof a particle number N , whose probability distributionis given by P ( N ). In experimental measurements, typi-cal detectors cannot count the particle number N in eachevent accurately due to the miss of particle observationand miscounting of ghost tracks. Therefore, the observedparticle number n by the detector in an event is generallydifferent from the true particle number N , and its event-by-event distribution denoted by ˜ P ( n ) is also modifiedfrom P ( N ). The goal of the detector-response correctionis to obtain the cumulants of P ( N ) from experimentally-observed distribution ˜ P ( n ) and the knowledge on the de-tector.To deal with this problem, we assume that the prob-ability to observe n particles in an event with the trueparticle number N depends only on N . We denote thisprobability as R ( n ; N ); because R ( n ; N ) is a probabil-ity, it satisfies (cid:80) n R ( n ; N ) = 1. We refer to R ( n ; N ) asthe response matrix. Using R ( n ; N ), ˜ P ( n ) and P ( N ) arerelated with each other as˜ P ( n ) = (cid:88) N R ( n ; N ) P ( N ) . (1)In this study we consider the detector-response correctionwith Eq. (1). Here, the response matrix R ( n ; N ) containsall information on the property of the detector relevantto this problem. We note that Eq. (1) can describe notonly the efficiency loss, but also the effects of the ghosttracks.Although the detector-response correction of a single-particle distribution is considered in Eq. (1), the correc-tion with multi-particle distributions are usually neces-sary in heavy-ion collisions [24]. In this study, however,we basically concentrate on the single-particle distribu-tion for simplicity. The generalization of the procedureto multi-particle distributions is straightforward as dis-cussed in Appendix A. B. Notation
We denote the m th-order moment of ˜ P ( n ) as (cid:104)(cid:104) n m (cid:105)(cid:105) = (cid:88) n n m ˜ P ( n ) , (2)while that of P ( N ) is expressed as (cid:104) N m (cid:105) = (cid:88) N N m P ( N ) . (3) The m th-order cumulants of these distributions are de-noted as (cid:104)(cid:104) n m (cid:105)(cid:105) c and (cid:104) N m (cid:105) c , respectively. C. Formal solution
Substituting Eq. (1) into Eq. (2) one finds that the m th-order moment of ˜ P ( n ) is given by (cid:104)(cid:104) n m (cid:105)(cid:105) = (cid:88) N P ( N ) (cid:88) n n m R ( n ; N )= (cid:88) N P ( N ) R m ( N ) , (4)where R m ( N ) is defined as R m ( N ) = (cid:88) n n m R ( n ; N ) . (5)Because R ( n ; N ) with fixed N represents a probability, R m ( N ) is understood as the moments of R ( n ; N ).We then suppose that R m ( N ) is expanded as R m ( N ) = (cid:88) j =0 r mj N j . (6)Substituting Eq. (6) into Eq. (4), one obtains (cid:104)(cid:104) n m (cid:105)(cid:105) = r m + (cid:88) j =1 r mj (cid:104) N j (cid:105) . (7)Equation (7) is expressed in the matrix form as (cid:104)(cid:104) n (cid:105)(cid:105)(cid:104)(cid:104) n (cid:105)(cid:105) ... = r r ... + R (cid:104) N (cid:105)(cid:104) N (cid:105) ... , (8)with the matrix R = ( r mj ). If R is a regular matrix, byapplying R − to Eq. (8) from left one obtains (cid:104) N (cid:105)(cid:104) N (cid:105) ... = R − (cid:104)(cid:104) n (cid:105)(cid:105)(cid:104)(cid:104) n (cid:105)(cid:105) ... − R − r r ... . (9)In Eq. (9), the moments of the true distribution (cid:104) N m (cid:105) are expressed by the experimentally-observed moments (cid:104)(cid:104) n m (cid:105)(cid:105) together with the parameters r mj characterizingthe property of the detector. Because the cumulants of P ( N ) are obtained from (cid:104) N m (cid:105) [1], the detector-responsecorrection of the cumulants is carried out with Eq. (9).In this argument, one can use alternative expansions of R m ( N ) instead of Eq. (6), which lead to the same finalresult but might lead to an enhancement of numericalstability. This strategy is discussed in Appendix B. D. Exactly solvable models
Although the above argument formally solves the prob-lem of the detector-response correction, Eq. (9) is notquite useful when the matrix R is not closed at some or-der. First of all, the inverse matrix of R is not determinedin general in this case. Second, even if the inverse ma-trix were obtained, Eq. (9) requires (cid:104)(cid:104) n m (cid:105)(cid:105) up to infinitelyhigher orders. In realistic situations, however, the mo-ments (cid:104)(cid:104) n m (cid:105)(cid:105) accessible with a reasonable statistics arelimited typically to m (cid:46) R ( n ; N ) satisfying this conditionare the cases that R m ( N ) is given by a m th-order poly-nomial, i.e. R m ( N ) = m (cid:88) j =0 r mj N j . (10)In this case, the matrix R in Eq. (8) has a lower-triangular form R = r · · · r r · · · r r r . . .... ... . . . . (11)The inverse of a lower-triangular matrix is obtainedorder-by-order, and R − is also lower triangular. Sub-stituting the lower-triangular form of R − into Eq. (9),one finds that (cid:104) N m (cid:105) depends only on (cid:104)(cid:104) n l (cid:105)(cid:105) for l ≤ m .The binomial model with R bin ( n ; N ) = B ( n ; p, N ) [1],with the binomial distribution B ( n ; p, N ) = p n (1 − p ) N − n N ! n !( N − n )! (12)corresponds to this case. In fact, all the cumulants of thebinomial distribution is proportional to N , (cid:104) n m (cid:105) c , binomial = ξ m ( p ) N, (13)where the coefficients ξ m ( p ) depend only on p [24]. Con-verting Eq. (13) into moments, one immediately findsthat R m ( N ) in this case is given by a m th-order poly-nomial as in Eq. (10). In this case, Eq. (9) reproducesthe formulas of the efficiency correction in the binomialmodel [1].Other examples satisfying Eq. (10) are the binomialmodel but the probability p is fluctuating event byevent [31], R G ( n ; N ) = (cid:90) dpG ( p ) B ( n ; p, N ) , (14)where G ( p ) is a probability distribution satisfying (cid:82) dpG ( p ) = 1. The moments R m ( N ) of Eq. (14) satisfyEq. (10) for arbitrary forms of G ( p ), as discussed in Ap-pendix C. The detector-response correction for R G ( n ; N )thus is handled with Eq. (9) exactly. The beta-binomialdistribution, which is obtained by G ( p ) = B ( p ; a, b ) withthe beta distribution B ( p ; a, b ), belongs to this case (seeAppendix D).Another interesting response matrix is the oneparametrized by the hypergeometric distribution as R HG ( n ; N ) = H ( n ; N, X, Y ) , (15) where the hypergeometric distribution H ( n ; N, X, Y ) isdefined in Appendix D. As shown in Appendix D, themoments of Eq. (15) is given in the form in Eq. (10).Therefore, the detector-response correction for Eq. (15)is also carried out exactly with Eq. (9).
E. Truncation
When the expansion of R m ( N ) is not closed, one mustintroduce an approximation to deal with the detector-response correction. A simple approximation is a trunca-tion of the expansion Eq. (6) at L th order, R m ( N ) = L (cid:88) j =0 r mj N j . (16)Using Eq. (16), one obtains a closed formula up to the L th order (cid:104)(cid:104) n (cid:105)(cid:105)(cid:104)(cid:104) n (cid:105)(cid:105) ... (cid:104)(cid:104) n L (cid:105)(cid:105) = r r ... r L + R (cid:104) N (cid:105)(cid:104) N (cid:105) ... (cid:104) N L (cid:105) , (17)where R is a L × L matrix. When R is a regular matrix,by applying R − from left Eq. (17) enables one to carryout the detector-response correction up to the L th orderusing the experimental data on (cid:104)(cid:104) n m (cid:105)(cid:105) for m ≤ L .Of course, this analysis can be justified only when thetruncated formula Eq. (16) well reproduces the functionalform of R m ( N ). When R ( n ; N ) is given by an analyticform, the effect of the truncation would be estimated an-alytically. When one considers the response matrices ofrealistic detectors, they are usually estimated by MonteCarlo simulations such as GEANT [32], which providethe moments R m ( N ) with statistical errors. In this case,one may perform fits to R m ( N ) with Eq. (16). The use ofEq. (17) would be justified as long as these fits reproduce R m ( N ) within statistics. The detector-response correc-tion of realistic detectors will be discussed in Sec. III inmore detail. III. PRACTICAL ANALYSIS
In this section, we discuss the detector-response correc-tion of realistic detectors whose response matrix R ( n ; N )are not given by an analytic form. In the following, weconsider the use of the approximation with the truncationdiscussed in Sec. II E.The form of R ( n ; N ) of realistic detectors is usually es-timated by Monte Carlo simulations such as GEANT [32].The simulations provide the moments R m ( N ) with statis-tical errors. In this case, the coefficients r mj in Eq. (16)are determined by the fits to R m ( N ) obtained by the sim-ulation. Using r mj thus obtained, the correction can becarried out with Eq. (17).Because we do not know the true distribution of N inrealistic situations, in the Monte Carlo simulations onemay assume a presumed “true” distribution P MC ( N ). Aproblem here is that the quality and result of the fits to R m ( N ) depend on the form of P MC ( N ) and the numberof the Monte Carlo events, N event . The validity of the fitswould be checked by setting N event to the same value asthe statistics of the experimental data. When the valueof chi-square, χ / ndf, of these fits are close to unity withthis statistics, there are no reasons to reject the use ofEq. (17). Next, the fitting results of r mj can also de-pend on the form of P MC ( N ). This suggests that onemust check the sensitivity of the fit results on the formof P MC ( N ), or perform an iterative procedure as follows:1. Generate R ( n ; N ) by a Monte-Carlo simulationwith a presumed distribution P MC ( N ).2. Perform fits to R m ( N ) with Eq. (16). One thenobtains r mj for m, j ≤ L . Together with the exper-imental results on (cid:104)(cid:104) n m (cid:105)(cid:105) , one obtains the correctedmoments (cid:104) N m (cid:105) .3. If (cid:104) N m (cid:105) thus obtained have large deviations fromthe moments of P MC ( N ), replace P MC ( N ) with theone consistent with (cid:104) N m (cid:105) obtained in the abovestep, and take the analysis from the top again.4. Repeat this iteration until P MC ( N ) is consistentwith (cid:104) N m (cid:105) obtained by the correction.It, however, is expected that the result of the fits are in-sensitive to P MC ( N ), especially on the cumulants higherthan the second order. The use of the Gaussian distri-bution with the mean and variance obtained by the cor-rection for P MC ( N ) would be sufficient for this analysis.It is also expected that a few iterations are enough forconvergence.Finally, we comment on the error analysis. First, inthe detector-response correction with Eq. (17), it is im-portant to reflect the correlation between the errors of r mj to the final result appropriately. An automatic wayto include the correlation is the use of the bootstrap orjackknife analysis with the successive generation of MonteCarlo events. Second, in the present method it is possibleto reduce the errors of r mj by increasing N event indepen-dently of the statistics of (cid:104)(cid:104) n m (cid:105)(cid:105) . In fact, in the nextsection we will see that the suppression of the error of r mj is effective in reducing the error of the final result.With increasing N event , however, the χ / ndf of the fits to R m ( N ) with Eq. (16) will eventually become unaccept-ably large. In this case, the analysis with the truncationloses its validity. In this sense, this analysis has an upperlimit of the resolution. Third, the effect of the truncationcan be estimated by comparing the corrected results atthe L and ( L + 1)th orders. Such analyses would requirelarge statistics, but are desirable for a proper estimate onthe systematic uncertainty of the analysis. IV. TEST ANALYSIS 1: EXACT MODELS
In this and next sections, we perform test analyses forthe detector-response correction discussed in Sec. II withtoy models for R ( n ; N ), and show that the correctionsare carried out successfully in these cases.In this section, we first perform test analyses for theresponse matrices which can be solved exactly discussedin Sec. II D. We consider two non-binomial models for FIG. 1. Correlation between n and N on the sample events,i.e. the magnitude of R ( n ; N ) P ( N ), for the response matrices R HG ( n ; N ) (hypergeometric) and R β ( n ; N ) (beta-binomial)with p = X/Y = 0 . Y = 140. R ( n ; N ) parametrized by the hypergeometric and beta-binomial distributions as R HG ( n ; N ) = H ( n ; N, X, Y ) , (18) R β ( n ; N ) = β ( n ; N, X, Y − X ) , (19)where the hypergeometric and beta-binomial distribu-tions, H ( n ; N, X, Y ) and β ( n ; N, a, b ), are defined in Ap-pendix D. The response matrices parametrized by thesedistributions are studied in Ref. [27] as examples thatthe binomial model fails in obtaining the true cumulants,and are good starting points for the check of the newmethod. Equations (18) and (19) approach the binomialmodel R bin ( n ; N ) = B ( n ; p, N ) in the Y → ∞ limit withfixed p = X/Y , while the distribution of n in R HG ( n ; N )( R β ( n ; N )) is narrower (wider) than the binomial dis-tribution with finite Y . As discussed in Appendix D,the values of r mj in Eq. (6) are obtained analytically for R HG ( n ; N ) and R β ( n ; N ).The procedure of the test analysis is as follows. Wefirst generate sample events of N by assuming the Poissondistribution for P ( N ) with (cid:104) N (cid:105) = 40. We then specifythe value of n for each sample event randomly accord-ing to the probability R HG ( n ; N ) or R β ( n ; N ). This al-lows one to obtain the moments (cid:104)(cid:104) n m (cid:105)(cid:105) . These momentsare used for the correction in Eq. (9). To proceed thecorrection, we take the following two different analyses.First, because the values of r mj are analytically knownfor R HG ( n ; N ) and R β ( n ; N ), we perform the correctionwith these values. Besides this analysis, as a second op-tion, we analyze (cid:104) N m (cid:105) with the values of r mj determinedby the fits to R m ( N ) obtained on the sample events withstatistical errors. The second analysis supposes the cor-rection of realistic detectors, of which the response matrixis obtained only stochastically.In Fig. 1, we show the correlation between n and N onthe 10 sample events by plotting the two-dimensionalhistogram as a function of n and N for the hypergeo-metric ( R HG ( n ; N )) and beta-binomial ( R β ( n ; N )) dis-tributions with p = 0 . Y = 140. (This plot thusrepresents the magnitude of R ( n ; N ) P ( N ), and is usu-ally called the “response matrix” in literature for sim-plicity.) One finds from the figure that the distributionsare clearly different between the two response matrices;the width of n with fixed N is narrower for R HG ( n ; N )than R β ( n ; N ).In Fig. 2, we show the cumulants of the response matrix FIG. 2. Cumulants of the response matrix C m ( N ) for R HG ( n ; N ) and R β ( n ; N ) obtained on 10 sample events with p = 0 . Y = 140. The dashed lines show the analytic val-ues, while the dotted lines represent the fitting results with m th-order polynomial. C m ( N ) defined by C ( N ) = R ( N ) , C ( N ) = R ( N ) − ( R ( N )) , (20)and so forth, for R HG ( n ; N ) and R β ( n ; N ) obtained on10 sample events with p = 0 . Y = 140 for m ≤
4. The dashed lines show the analytic values, while thedotted lines are the fitting results with the m th-orderpolynomial. From these fits one obtains the values of r mj .In Fig. 3, we show the corrected values of the cumu-lants (cid:104) N m (cid:105) c for m ≤ p = 0 . Y . The left (right) panel shows the results for R HG ( n ; N )( R β ( n ; N )). The triangles represent the results obtainedwith the analytic values of r mj , while the results obtainedwith r mj determined by the fits to R m ( N ) are shown bysquares. 10 sample events are used to obtain (cid:104)(cid:104) n m (cid:105)(cid:105) inboth analyses, while r mj in the latter analysis are ob-tained with 10 sample events. Errors are estimated byrepeating the same simulation 100 times. One finds fromthe figure that the corrected cumulants (cid:104) N m (cid:105) c are con-sistent with the true value, (cid:104) N m (cid:105) c = 40 shown by thedashed line, within statistics for all values of Y in bothanalyses. In Fig. 3, the uncorrected cumulants, (cid:104)(cid:104) n m (cid:105)(cid:105) c ,are shown by filled circles. We also show the results ofthe efficiency correction with the binomial model with p = 0 . r mj are determined by the fits, althoughthe statistics to determine r mj is one order larger thanthat for (cid:104)(cid:104) n m (cid:105)(cid:105) . This suggests that the suppression of theuncertainty of r mj is crucial in reducing the error of thefinal results. FIG. 3. Cumulants obtained by the detector-response correc-tion, (cid:104) N m (cid:105) , up to the fourth order with p = 0 . Y for R HG ( n ; N ) (left) and R β ( n ; N ) (right). The resultsobtained with the analytic (fitted) values of r mj are shownby triangles (squares). The corrected values agree with thetrue cumulants (cid:104) N m (cid:105) c = 40 shown by the dashed line withinstatistics. The uncorrected cumulants (cid:104)(cid:104) n m (cid:105)(cid:105) c and the cor-rected results in the binomial model are also shown by circlesand stars, respectively. Finally, we note that the fitting results of C m ( N ) inFig. 2 have significant deviations from the analytic val-ues for N (cid:38)
60. Nevertheless, the final results obtainedwith these fits reproduce the true values within statistics.This result shows that the detector-response correction iscarried out appropriately even if the fits do not reproduce R m ( N ) in the range of N at which P ( N ) is small. V. TEST ANALYSIS 2:MULTIPLICITY-DEPENDENT EFFICIENCY
Next, we perform a test analysis of the detector-response correction for the response matrix which cannotbe solved exactly. As such an example, we consider theresponse of a detector having a multiplicity-dependent ef-ficiency. We consider the binomial distribution but theefficiency is dependent on N , i.e. R MD ( n ; N ) = B ( n ; p ( N ) , N ) . (21) FIG. 4. Cumulants C m ( N ) of the response matrix R MD ( n ; N ) with p = 0 . ε = 0 .
002 up to the fourthorder obtained on 10 sample events. The dashed lines showthe results of the fits by the fourth order polynomial. In typical detectors, the efficiency decreases with increas-ing multiplicity ( N ) [9, 13]. To model this behavior weassume p ( N ) = 0 . − ε ( N − (cid:104) N (cid:105) ) , (22)with ε > One can analytically show that the m th-order moment R m ( N ) of R MD ( n ; N ) with Eq. (22) is given by the 2 m th-order polynomial. Therefore, the detector-response cor-rection in this case cannot be solved exactly by the pro-cedure in Sec. II D. In the following, we use the truncatedformulas discussed in Sec. II E with L = 4.For the test analysis, we generate the sample eventsof N assuming the Poisson distribution with (cid:104) N (cid:105) = 40for P ( N ). The value of n in each sample event is thenspecified according to R MD ( n ; N ). 10 sample events aregenerated in this way. Figure 4 shows the cumulants ofthe response matrix, C m ( N ), for m ≤ p = 0 . ε = 0 . R m ( N ) with the fourth-order polynomialEq. (16). The results of the fits in terms of C m ( N ) areshown by the dashed lines in Fig. 4 . As shown in thefigure, we have χ / ndf (cid:39) The distribution of n of realistic detectors with fixed N would bewider or narrower than the binomial distribution. To model sucha behavior with the multiplicity-dependent efficiency Eq. (22),one may, for example, employ the response matrix R ( n ; N ) = β ( n ; N, p ( N ) Y, Y − p ( N ) Y ) or R ( n ; N ) = H ( n ; N, p ( N ) Y, Y ). We also tested the fits to C m ( N ), instead of R m ( N ), by Eq. (16).We have checked that the results converted to R m ( N ) agrees withthose obtained by the fit to R m ( N ) directly to a good accuracy. FIG. 5. Filled circles show the corrected cumulants (cid:104) N m (cid:105) c obtained by the truncation analysis for R MD ( n ; N ) as func-tions of ε . The corrected cumulants reproduce the true value (cid:104) N m (cid:105) c = 40 for all ε within statistics. The open circles showthe corrected results in the binomial model. The results of the corrected cumulants, (cid:104) N m (cid:105) c , areshown in Fig. 5 by filled circles as a function of ε for m ≤
4. Errors are estimated by the same procedureas in the previous section. The figure shows that thecorrected results reproduce (cid:104) N m (cid:105) c = 40 shown by thedashed line in each panel within statistics for all values of ε . In Fig. 5, we also show the cumulants obtained by theefficiency correction in the binomial model by the opencircles. These results fail in reproducing the true cumu-lants even at m = 1. Although the results in the binomialmodel are consistent with the true value for m = 4, thisagreement would be accidental. We note that the typicalvalue of ε for protons at the STAR detector is ε (cid:39) . √ s NN = 7 . ε (cid:39) . √ s NN = 200 GeVin Au+Au collisions [33]. The maximal value of ε inFig. 5 is about one order larger than this value.As discussed in Sec. III, in realistic experimental anal-ysis we do not know the true distribution P ( N ). Thismeans that the distribution P MC ( N ) used in the MonteCarlo simulation to determine R ( n ; N ) is in general dif-ferent from P ( N ). In order to test the detector-responsecorrection in such a case, we performed one more testanalysis as follows: We use the Poisson distribution with (cid:104) N (cid:105) = 40 for P ( N ) to determine (cid:104)(cid:104) n m (cid:105)(cid:105) as in the aboveanalysis, but the values of r mj are determined by thefits to R m ( N ) generated by the Gauss distribution for P MC ( N ) with (cid:104) N (cid:105) = (cid:104) N (cid:105) c = 40. We checked that thefinal results in this analysis agrees with those shown inFig. 5 within statistics even in this case. This result sug- We note that these numbers are the slope with respect to themultiplicity used for the centrality definition, which are takeninto account in the experimental analysis at STAR by using theCentraltiy Bin Width Correction [2]. There is no study thatdiscusses the slope with respect to the net-particle itself, yet. gests that the detector-response correction with the strat-egy in Sec. III is applicable to realistic analyses.
VI. SUMMARY
In this paper, we proposed a new procedure to carryout the detector-response correction (including efficiencycorrection) of the higher order cumulants in the event-by-event analysis. This method provides exact formu-las for various non-binomial response matrices, includ-ing Eqs. (18) and (19) parametrized by the hypergeo-metric and beta-binomial distributions. Introducing anapproximation with the truncation, this method can dealwith the correction of realistic detectors. The correctionin this method is demonstrated explicitly for three non-binomial response matrices. We showed that the truecumulants are obtained within statistics not only for theexactly-solvable cases but also a case that the approxi-mation with the truncation is necessary. Although weconcentrated on the correction for the single-variable dis-tribution throughout this manuscript, the extension todeal with the multi-variable distribution is also possibleas discussed in Appendix A.” “We thank M. Gazdzicki, R. Holzmann, A. Rustamov,M. Szala, and N. Xu for constructive discussions. M. K.thanks K. Redlich and B. Friman for inviting him tothe workshop “Constraining the QCD Phase Boundarywith Data from Heavy Ion Collisions” (Feb. 12-14, 2018,Darmstadt, Germany) and fruitful discussions during theworkshop.
Appendix A: Detector-response correction formulti-variable distribution
In this appendix, we discuss the extension of thedetector-response correction to multi-variable distribu-tion functions. This extension is necessary for the analy-sis of the net-particle number cumulants, because the netnumber is given by the difference of particle and anti-particle numbers. Also, the effects of the momentum andazimuthal angle dependence of the detector’s responsecan be described by the multi-variable distribution [24].For a simple illustration, we consider a distributionfunction of two particle species, P ( N , N ), where N and N are the numbers of the particles in each event. Wethen denote the numbers of observed particles in eachevent as n and n , respectively, and its event-by-eventdistribution as ˜ P ( n , n ).Similar to Eq. (1), we assume that these distributionfunctions are related with each other as˜ P ( n , n ) = (cid:88) N ,N R ( n , n ; N , N ) P ( N , N ) , (A1)where the response matrix R ( n , n ; N , N ) satisfies thenormalization condition (cid:80) n ,n R ( n , n ; N , N ) = 1.Next, similar to Eq. (5) we consider the moments of R ( n , n ; N , N ). In this case, because R ( n , n ; N , N )has two variables, we must consider the mixed moments, R m m ( N , N ) = (cid:88) n ,n n m n m R ( n , n ; N , N ) . (A2) We then Taylor expand R m m ( N , N ), R m m ( N , N ) = (cid:88) j ,j r m m ; j j N j N j . (A3)Using Eq. (A3), the mixed moments of observed particlenumbers (cid:104)(cid:104) n m n m (cid:105)(cid:105) = (cid:88) n ,n n m n m ˜ P ( n , n ) (A4)are given by (cid:104)(cid:104) n m n m (cid:105)(cid:105) = r m m ;00 + (cid:88) j ,j r m m ; j j (cid:104) N j N j (cid:105) . (A5)We note that Eqs. (A3) and (A4) corresponds to Eqs. (6)and (7), respectively. In the matrix form, Eq. (A5) isexpressed as (cid:104)(cid:104) n (cid:105)(cid:105)(cid:104)(cid:104) n (cid:105)(cid:105)(cid:104)(cid:104) n (cid:105)(cid:105)(cid:104)(cid:104) n n (cid:105)(cid:105)(cid:104)(cid:104) n (cid:105)(cid:105) ... = r r r r r ... + R mult (cid:104) N (cid:105)(cid:104) N (cid:105)(cid:104) N (cid:105)(cid:104) N N (cid:105)(cid:104) N (cid:105) ... , (A6)where the matrix R mult is composed of r m m ; j j . Byinversely solving Eq. (A6) as in Eq. (9), one obtains theformulas to obtain the (mixed-)moments of the true dis-tribution function (cid:104) N m N m (cid:105) . The cumulants are thenconstructed from these moments; see Refs. [24, 34–36],for the construction of the cumulants with multi-particlespecies. When the expansion Eq. (A3) is not closed, onemay truncate the expansion as in Eq. (16); the truncationat j + j ≤ L allows one to carry out the correction upto L th order.This procedure can be generalized to the case withmore than two-particle species in a straightforward man-ner. Appendix B: Alternative expansions of R m ( N ) In this appendix, we consider modification of the pro-cedure in Sec. II with the use of alternative expansions of R m ( N ). These procedures give the same final result inprinciple, but might be effective in suppressing the accu-mulation of numerical errors in practical analyses.First, instead of Eq. (16) we consider the Taylor ex-pansion R m ( N ) = L (cid:88) j =0 ¯ r mj ( N − N ) j (B1)at N = N with an arbitrary number N . SubstitutingEq. (B1) into Eq. (4), one obtains (cid:104)(cid:104) n (cid:105)(cid:105)(cid:104)(cid:104) n (cid:105)(cid:105) ... (cid:104)(cid:104) n L (cid:105)(cid:105) = ¯ r ¯ r ...¯ r L + ¯ R (cid:104) N − N (cid:105)(cid:104) ( N − N ) (cid:105) ... (cid:104) ( N − N ) L (cid:105) , (B2)with ¯ R = (¯ r mj ). By choosing N = (cid:104) N (cid:105) , Eq. (B2) allowsone to obtain the central moments (cid:104) ( N −(cid:104) N (cid:105) ) m (cid:105) directly.We stress that Eq. (B1) represents the same functionas Eq. (16) with an appropriate replacement between r mj and ¯ r mj . Therefore, if these parameters are determinedaccurately Eqs. (17) and (B2) give the same final result.However, they can give different results within numericalprecision in practice. The use of Eq. (B2) would be ad-vantageous in reducing the numerical error, as the centralmoments are more closely related to cumulants than thestandard moments.Second, it is also possible to use expansions motivatedby the factorial moments. One may expand R m ( N ) as R m ( N ) = f m + f m N + f m N ( N −
1) + · · · + f mj N ( N − · · · ( N − j + 1) + · · · . (B3)Substituting this expansion into Eq. (4), one obtains (cid:104)(cid:104) n (cid:105)(cid:105)(cid:104)(cid:104) n (cid:105)(cid:105) ... = f f ... + F (cid:104) N (cid:105) f (cid:104) N (cid:105) f ... , (B4)where F = ( f mj ) and (cid:104) N j (cid:105) f = (cid:104) N ( N − · · · ( N − j + 1) (cid:105) , (B5)are the factorial moments of P ( n ). Equation (B4) pro-vides formulas to obtain factorial moments (cid:104) N m (cid:105) f di-rectly.One can also expand the factorial moments of R ( n ; N )as R f ,m ( N ) = (cid:88) n n ( n − · · · ( n − m + 1) R ( n ; N )= ¯ f m + ¯ f m N + ¯ f m N ( N −
1) + · · · + ¯ f m N ( N − · · · ( N − j + 1) + · · · , (B6)which enables us to connect (cid:104) N m (cid:105) f with the factorial mo-ments of ˜ P ( N ), (cid:104)(cid:104) n m (cid:105)(cid:105) f , directly as (cid:104)(cid:104) n (cid:105)(cid:105) f (cid:104)(cid:104) n (cid:105)(cid:105) f ... = ¯ f ¯ f ... + ¯ F (cid:104) N (cid:105) f (cid:104) N (cid:105) f ... , (B7)with ¯ F = ( ¯ f mj ). In the binomial model, ¯ F is given by adiagonal matrix [34]. When R ( n ; N ) is well approximatedby the binomial model, therefore, the numerical analysisof the inverse matrix of ¯ F would be more stable than thatof ¯ R . Appendix C: Binomial distribution with fluctuatingprobability
In this appendix, we calculate the moments of the bi-nomial distribution Eq. (12) but the probability p is fluc-tuating, i.e. the probability distribution given by P N ( n ) = (cid:90) dpG ( p ) B ( n ; p, N ) , (C1) We thank R. Holzmann for pointing out this property of thematrix ¯ F . where G ( p ) satisfies (cid:82) dpG ( p ) = 1.To obtain the moments of P N ( n ), (cid:104) n m (cid:105) = (cid:88) n n m P N ( n ) , (C2)it is convenient to first calculate their factorial moments,which are given by (cid:104) n m (cid:105) f = (cid:88) n n ( n − · · · ( n − m + 1) P N ( n )= (cid:90) dpG ( p ) (cid:88) n n ( n − · · · ( n − m + 1) B ( n ; p, N )= (cid:90) dpG ( p ) N ( N − · · · ( N − m + 1) p m = N ( N − · · · ( N − m + 1) (cid:104) p m (cid:105) G , (C3)where (cid:104) p m (cid:105) G = (cid:82) dpp m G ( p ) is the moments of G ( p ) andin the last equality we used the relation (cid:88) n n ( n − · · · ( n − m + 1) B ( n ; p, N )= N ( N − · · · ( N − m + 1) p m . (C4)The moments of P N ( n ) is then obtained by convertingEq. (C3) [34]. Explicit results up to the fourth order are (cid:104) n (cid:105) = (cid:104) n (cid:105) f = (cid:104) p (cid:105) G N, (C5) (cid:104) n (cid:105) = (cid:104) n (cid:105) f + (cid:104) n (cid:105) f = (cid:104) p − p (cid:105) G N + (cid:104) p (cid:105) G N , (C6) (cid:104) n (cid:105) = (cid:104) n (cid:105) f + 3 (cid:104) n (cid:105) f + (cid:104) n (cid:105) f = (cid:104) p − p + 2 p (cid:105) G N + 3 (cid:104) p − p (cid:105) G N + (cid:104) p (cid:105) G N , (C7) (cid:104) n (cid:105) = (cid:104) n (cid:105) f + 7 (cid:104) n (cid:105) f + 6 (cid:104) n (cid:105) f + (cid:104) n (cid:105) f = (cid:104) p − p + 12 p − p (cid:105) G N + (cid:104) p − p + 11 p (cid:105) G N + 6 (cid:104) p − p (cid:105) G N + (cid:104) p (cid:105) G N . (C8)The same manipulation can be repeated for arbitraryhigher orders. From this derivation, it is clear that the m th-order moments of P N ( n ) are given by the m th-orderpolynomial of N . Appendix D: Hypergeometric and beta-binomialdistributions
In this appendix we summarize the definitions andproperties of the hypergeometric and beta-binomial dis-tributions.We define the hypergeometric and beta-binomial dis-tributions, H ( n ; N, X, Y ) and β ( n ; N, a, b ), as H ( n ; N, X, Y )= X ! n !( X − n )! ( Y − X )!( N − n )!( X − n )! N !( Y − N )! Y ! , (D1)and β ( n ; N, a, b ) = (cid:90) dp B ( p ; a, b ) B ( n ; p, N ) , (D2)with the beta distribution B ( p ; a, b ) = p a (1 − p ) b /B ( a, b ) , (D3)where B ( a, b ) is the beta function required for the nor-malization (cid:82) dp B ( p ; a, b ) = 1.The hypergeometric and beta distributions are given inurn models as follows. First, we consider N w white ballsand N b black balls in an urn, and draw N balls fromthe urn without returning balls to the urn in each draw.Then, the probability distribution of the number of whiteballs, n , is given by the hypergeometric distribution as H ( n ; N, N w , N tot ) with N tot = N w + N b . Next, we con-sider successive draws of the balls from the urn. In eachdraw, when one draws a white (black) ball, two white (black) balls are returned to the urn. After repeatingthis procedure N times, the probability distribution todraw n white balls in total is given by the beta-binomialdistribution as β ( n ; N, N w , N b ). The distributions of n in both urn models become close to the binomial dis-tribution B ( n ; p, N ) in the limit N tot → ∞ with fixed p = N w /N tot , where p represents the probability to drawa white ball in a draw. This means that H ( n ; N, X, Y )and β ( n ; N, a, b ) approach the binomial distribution inthe limit Y → ∞ with fixed X/Y and a → ∞ with fixed a/ ( a + b ), respectively.The cumulants of Eqs. (D1) and (D2), (cid:104) n m (cid:105) HGc and (cid:104) n m (cid:105) β c , respectively, up to the fourth order are given by (cid:104) n (cid:105) HGc = N p, (D4) (cid:104) n (cid:105) HGc = N ( N − Y )( − p ) p − Y , (D5) (cid:104) n (cid:105) HGc = N (2 N − N Y + Y )( − p ) p ( − p )( − Y )( − Y ) , (D6) (cid:104) n (cid:105) HGc = [ N ( N − Y )( − p ) p (6 N ( − − − p ) p + Y (1 + 5( − p ) p )) − N Y ( − − − p ) p + Y (1 + 5( − p ) p ))+( − Y ) Y (1 + Y (1 + 6( − p ) p )))] / (cid:2) ( − Y )( − Y )( − Y ) (cid:3) , (D7) (cid:104) n (cid:105) β c = N p, (D8) (cid:104) n (cid:105) β c = − N ( N + a + b )( − p ) p a + b , (D9) (cid:104) n (cid:105) β c = N ( N + a + b )(2 N + a + b )( − p ) p ( − p )(1 + a + b )(2 + a + b ) , (D10) (cid:104) n (cid:105) β c = − (( N ( N + a + b )( − p ) p (6 N (1 + 6( − p ) p + ( a + b )(1 + 5( − p ) p ))+6 N ( a + b )(1 + 6( − p ) p + ( a + b )(1 + 5( − p ) p )) + ( a + b )(1 + a + b )( − a + b )(1 + 6( − p ) p )))) / ((1 + a + b ) (2 + a + b )(3 + a + b ))) . (D11)As one finds from Eq. (D5) (Eq. (D9)), the hypergeomet-ric (beta-binomial) distribution is narrower (wider) thanthe binomial distribution with (cid:104) n (cid:105) c = p (1 − p ) N .As the m th-order cumulants are given by the m th-order polynomial in Eqs. (D4)–(D11), the m th-order momentsof Eqs. (D1) and (D2) are also given by m th-order poly-nomial. The values of r mj are obtained as the coefficientsof the moments. [1] M. Asakawa and M. Kitazawa, Prog. Part. Nucl. Phys. , 299 (2016), arXiv:1512.05038 [nucl-th].[2] X. Luo and N. Xu, Nucl. Sci. Tech. , 112 (2017),arXiv:1701.02105 [nucl-ex].[3] S. Ejiri, F. Karsch, and K. Redlich, Phys. Lett. B633 ,275 (2006), arXiv:hep-ph/0509051 [hep-ph].[4] M. A. Stephanov, Phys. Rev. Lett. , 032301 (2009),arXiv:0809.3450 [hep-ph].[5] M. Asakawa, S. Ejiri, and M. Kitazawa, Phys. Rev. Lett. , 262301 (2009), arXiv:0904.2089 [nucl-th].[6] L. Adamczyk et al. (STAR), Phys. Rev. Lett. , 032302(2014), arXiv:1309.5681 [nucl-ex].[7] L. Adamczyk et al. (STAR), Phys. Rev. Lett. , 092301(2014), arXiv:1402.1558 [nucl-ex].[8] A. Adare et al. (PHENIX), Phys. Rev.
C93 , 011901(2016), arXiv:1506.07834 [nucl-ex]. [9] L. Adamczyk et al. (STAR), (2017), arXiv:1709.00773[nucl-ex].[10] B. Abelev et al. (ALICE), Phys. Rev. Lett. , 152301(2013), arXiv:1207.6068 [nucl-ex].[11] A. Rustamov (ALICE),
Proceedings, 26th InternationalConference on Ultra-relativistic Nucleus-Nucleus Colli-sions (Quark Matter 2017): Chicago, Illinois, USA,February 5-11, 2017 , Nucl. Phys.
A967 , 453 (2017),arXiv:1704.05329 [nucl-ex].[12] M. Gazdzicki, talk at “Constraining the QCD PhaseBoundary with Data from Heavy Ion Collisions” (Darm-stadt, Germany, Feb. 12-14, 2018).[13] R. Holzmann (HADES), Talk at “26th International Con-ference on Ultra-relativistic Nucleus-Nucleus Collisions(Quark Matter 2017)” (Chicago, Illinois, USA, Feb. 5-11, 2017). [14] H. Sako et al. , “White paper for a Future J-PARC Heavy-Ion Program (J-PARC-HI),” http://asrc.jaea.go.jp/soshiki/gr/hadron/jparc-hi/ .[15] R. Rapp et al. , Lect. Notes Phys. , 335 (2011).[16] D. Blaschke et al. , Eur. Phys. J. A , 299 (2016).[17] M. Kitazawa and M. Asakawa, Phys. Rev. C86 ,024904 (2012), [Erratum: Phys. Rev.C86,069902(2012)],arXiv:1205.3292 [nucl-th].[18] A. Bialas and R. B. Peschanski, Nucl. Phys.
B273 , 703(1986).[19] A. Bzdak and V. Koch, Phys. Rev.
C86 , 044904 (2012),arXiv:1206.4286 [nucl-th].[20] A. Bzdak and V. Koch, Phys. Rev.
C91 , 027901 (2015),arXiv:1312.4574 [nucl-th].[21] X. Luo, Phys. Rev.
C91 , 034907 (2015), arXiv:1410.3914[physics.data-an].[22] M. Kitazawa, Phys. Rev.
C93 , 044911 (2016),arXiv:1602.01234 [nucl-th].[23] T. Nonaka, T. Sugiura, S. Esumi, H. Masui, and X. Luo,Phys. Rev.
C94 , 034909 (2016), arXiv:1604.06212 [nucl-th].[24] T. Nonaka, M. Kitazawa, and S. Esumi, Phys. Rev.
C95 ,064912 (2017), arXiv:1702.07106 [physics.data-an].[25] R. Esha (STAR),
Proceedings, 26th International Con-ference on Ultra-relativistic Nucleus-Nucleus Collisions(Quark Matter 2017): Chicago, Illinois, USA, February5-11, 2017 , Nucl. Phys.
A967 , 457 (2017).[26] X. Luo (STAR),
Proceedings, 9th International Workshop on Critical Point and Onset of Deconfinement (CPOD2014): Bielefeld, Germany, November 17-21, 2014 , PoS
CPOD2014 , 019 (2015), arXiv:1503.02558 [nucl-ex].[27] A. Bzdak, R. Holzmann, and V. Koch, Phys. Rev.
C94 ,064907 (2016), arXiv:1603.09057 [nucl-th].[28] G. D’Agostini, Nucl. Instrum. Meth.
A362 , 487 (1995).[29] T. Adye, in
Proceedings, PHYSTAT 2011 Workshop onStatistical Issues Related to Discovery Claims in SearchExperiments and Unfolding, CERN,Geneva, Switzerland17-20 January 2011 , CERN (CERN, Geneva, 2011) pp.313–318, arXiv:1105.1160 [physics.data-an].[30] P. Garg, D. K. Mishra, P. K. Netrakanti, A. K. Mo-hanty, and B. Mohanty, J. Phys.
G40 , 055103 (2013),arXiv:1211.2074 [nucl-ex].[31] S. He and X. Luo, (2018), arXiv:1802.02911[physics.data-an].[32] S. Agostinelli et al. (GEANT4), Nucl. Instrum. Meth.
A506 , 250 (2003).[33] T. Nonaka, talk at “XII Workshop on Particle Cor-relations and Femtoscopy ” (Nikhef, Amsterdam, TheNetherlands, Jun. 12-16, 2017).[34] M. Kitazawa and X. Luo, Phys. Rev.
C96 , 024910 (2017),arXiv:1704.04909 [nucl-th].[35] E. A. De Wolf, I. M. Dremin, and W. Kittel, Phys. Rept. , 1 (1996), arXiv:hep-ph/9508325 [hep-ph].[36] W. Broniowski and A. Olszewski, Phys. Rev.