Efficiency correction for cumulants of multiplicity distributions based on track-by-track efficiency
Xiaofeng Luo Toshihiro Nonaka


 


xfluo@mail.ccnu.edu.cn
Key Laboratory of Quark & Lepton Physics (MOE) and Institute of Particle Physics, Central China Normal University, Wuhan 430079, China

tnonaka@rcf.rhic.bnl.gov
Key Laboratory of Quark & Lepton Physics (MOE) and Institute of Particle Physics, Central China Normal University, Wuhan 430079, China





We propose a simplified procedure for the experimental application of 
the efficiency correction on higher order cumulants in heavy-ion collisions. By using the track-by-track efficiency, we can eliminate possible bias arising from the average efficiencies calculated within the arbitrary binning 
of the phase space. Furthermore, the corrected particle spectra is no longer necessary for the average efficiency estimation and 
the time cost for the calculation of bootstrap statistical error can be significantly reduced. 

Efficiency correction for cumulants of multiplicity distributions based on track-by-track efficiency
    Toshihiro Nonaka
    December 30, 2023
====================================================================================================









§ INTRODUCTION


Higher order cumulants of event-by-event conserved charge fluctuations 
are important observables to search for the 
QCD critical point <cit.> in heavy-ion collisions.
In the beam energy scan (BES) program at RHIC, 
the STAR experiment has measured the cumulants up to the 4th order 
of the net-charge and net-kaon multiplicity 
distributions <cit.>,
and up to the 6th order of the net-proton multiplicity 
distribution <cit.>.
In the recent results of the net-proton 
multiplicity distribution, the non-monotonic behaviour 
of the fourth order fluctuation has been observed with 
respect to the collision energy <cit.>, 
which is quite similar to the theoretical prediction 
with the critical point <cit.>.  
Since there are still large uncertainties in low collision energies, 
the second phase of the Beam Energy Scan program (BES-II) will be carried out in 2019-2021 focusing 
on the collision energy of √(s_ NN)=7.7–19.6 GeV. 

One of the difficulties of measuring the higher order cumulants is the efficiency correction.
It is known that the value of cumulants are artificially changed from the true ones <cit.>,
due to the fact that detectors miss some particles with the probability called efficiency. 
Some analytical formulas have been proposed to correct the measured cumulants 
with the assumption that the response function of the efficiency follows the binomial distribution 
 <cit.>.
Recently, a few attempts to understand and correct for the possible non-binomial 
efficiencies have been discussed <cit.>, 
but there are still large systematic uncertainties arising from 
how to determine the detector-response functions.
The efficiency correction with the binomial assumption is thus still important. 
Hereafter, let us call it "binomial correction" for simplicity. In the binomial correction, particles are counted separately 
in the "efficiency bin" where the efficiency changes, 
which are substituted into the correction formulas 
with the corresponding value of efficiencies.
Due to the huge calculation cost with large number of efficiency bins
by using the correction formulas based on 
the factorial moments <cit.>, 
more efficient formulas have been proposed in which factorial 
cumulants are used in the derivation <cit.>.
Experimentally, single particle efficiencies can be computed 
by the MC detector simulations in terms of 
the various experimental observables like 
centrality, multiplicity, vertex position, 
transverse momentum, rapidity, azimuthal angle and so on.
Efficiency bins are then defined with respect to 
the track-wise variables in which the particles are counted. 
This binning is still arbitrary, which depends not only on the computing power 
for the detector simulations but also on 
the homogeneousness of the detector in the acceptance.


Therefore, additional systematic studies will be necessary 
on how many efficiency bins are enough. 
The most crucial thing is that we need to calculate the averaged 
efficiency at each efficiency bin with weighted by the true spectra. 
This indicates that the traditional efficiency correction based on average efficiency, can not be performed 
until the corrected spectra of identified particle is available. Fortunately, we found those difficulties can be overcome by using the so called track-by-track efficiency correction methods. 
By doing this, the corrected particle spectra is no longer needed. We can also reduce the potential systematic bias by using the average efficiency correction method and the calculation cost for bootstrap statistical errors calculations.











This paper is organized as follows. 
In Sec <ref>, the efficiency correction with the binomial model 
is introduced. The correction formulas based on factorial cumulants will be shown as well.
In Sec. <ref>, we discuss three difficulties 
of the efficiency correction in the experimental applications. 
We also use a numerical study to demonstrate one of the issues.
In Sec. <ref>, a simple solution using the 
track-by-track efficiency is shown, followed by the toy model 
to check the validity of the new method.









§ EFFICIENCY CORRECTION




 §.§ Cumulants and factorial cumulants

The m-th order cumulant of the probability distribution function P(N) 
is defined as 

    ⟨ N^m⟩_ c   =   ∂^m/∂θ^mK(θ)|_θ=0, 
     
    	K(θ)    =   ln∑_Ne^NθP(N) = ln⟨ e^Nθ⟩,

where K(θ) represents the cumulant generating function. 
Another quantity, factorial cumulants ⟨ N^m⟩_ fc 
are also defined as

    ⟨ N^m⟩_ fc   =   ∂^m/∂^mK_ f(s)|_s=1, 
    
    	K_ f(s)    =   ln⟨ s^N⟩,

with K_ f(s) being the factorial-cumulant generating function.
Cumulants and factorial cumulants are connected to each other, 
e.g, factorial cumulants are expressed in terms of cumulants: 

    ⟨ N^m⟩_ fc = ⟨ N(N-1)⋯(N-m+1)⟩_ c.





 §.§ Binomial model

Let us assume that the probability distribution function P(N) is 
observed as P̃(n) through detectors. 
Some of N generated particles are missed by the detectors, 
which leads to the observation of n particles (n≤ N).
This finite probability to measure particles 
characterized by the detector is called efficiency.
When the efficiencies for generated particles are independent each other, 
the detection process can be described 
by the binomial distribution B_ε,N(n):

    P̃(n)    =   ∑_NP(N)B_ε,N(n), 
     
    	B_ε,N(n)    =   N!/n!(N-n)!ε^n(1-ε)^N-n,
  
where ε represents the efficiency.
It is known that in this situation the relationship between
factorial cumulants of P(N) and P̃(n) 
is given by <cit.>

    ⟨ n^m⟩_ fc 
    	= ε^m⟨ N^m⟩_ fc.

Using Eq. (<ref>) and extending to the M multi-variable case 
P( N)=P(N_1,N_2,...,N_M),
the formulas up to the fourth order cumulant are shown below:


    <Q>_ c   =   q_(1,1)_ c, 
            
    
    <Q^2>_ c   =   q_(1,1)^2_ c 
    	+ q_(2,1)_ c - q_(2,2)_ c, 
            
    
    <Q^3>_ c   =   q_(1,1)^3_ c 
            + 3q_(1,1)q_(2,1)_ c - 3q_(1,1)q_(2,2)_ c
            + q_(3,1)_ c - 3q_(3,2)_ c
            + 2q_(3,3)_ c, 
    
    <Q^4>_ c   =   q_(1,1)^4_ c 
            + 6q_(1,1)^2q_(2,1)_ c - 6q_(1,1)^2q_(2,2)_ c 
    	+ 4q_(1,1)q_(3,1)_ c + 3q_(2,1)^2_ c
           + 3q_(2,2)^2_ c - 12q_(1,1)q_(3,2)_ c 
    	+ 8q_(1,1)q_(3,3)_ c -  6q_(2,1)q_(2,2)_ c
           + q_(4,1)_ c - 7q_(4,2)_ c 
            + 12q_(4,3)_ c - 6q_(4,4)_ c,

with q_(r,s) defined as

    q_(r,s) = 
            ∑_i=1^M (a_i^r/ε_i^s) n_i,


where M represents the number of efficiency bins, 
n_i represents the number of particles, 
and a_i represents the electric charge 
of particles in ith efficiency bin.
The concept of the efficiency bin started to be taken into account 
in Ref. <cit.> inspired by experimental requirement, 
which will be discussed in the next section. 



§ DIFFICULTIES IN THE EFFICIENCY CORRECTION


In this section, we clarify the three difficulties in the 
experimental application of the efficiency correction. 
First, the details of the experimental procedures 
of the efficiency correction is explained to point 
out the difficulties. 
Second, a simple toy model is used to demonstrate one of those. 



 §.§ Experimental application

Experimentally, single-particle efficiency can be determined
by the Monte-Carlo approach of the detector simulation 
with respect to track-wise variable like transverse momentum, 
rapidity and azimuthal angle. This kind of efficiency map can be computed precisely 
as long as the computing source permits.
Usually,  the efficiency map is divided into various efficiency bins based on the track-wise variables
and the averaged efficiency in each efficiency bin needs to be estimated by using the corrected spectra.
In the current analysis of the net-proton fluctuations 
at the STAR experiment, the proton identification method 
is different between low and high p_T regions, 
which leads to the step-like dependence of the efficiency 
as a function of p_T <cit.>.
In this case, particles need to be counted separately 
for the two p_T region in which the values of the efficiency 
are different, then the true cumulants 
in entire p_T regions can be reconstructed based on factorial moments <cit.>.

If the number of efficiency bins M is large, it is more efficient to apply the efficiency correction according to 
Eqs. (<ref>)–(<ref>), which is based on factorial cumulants. 
For the statistical errors estimation, 
the bootstrap method would be the realistic way. 
We can obtain a new distribution by random sampling 
from the original distribution with the same number of events to calculate cumulants.
This procedure is repeated with 100 times, 
and the standard deviation of the 100 cumulant values calculated from the new distributions
are taken as the statistical error. 
In order to take into account the correlation between different 
efficiency bins, the sampling would be performed based on 
the M dimensional histogram.

Below are three main difficulties in the experimental implementation  
of the efficiency correction.

	
  * We expected that it can
	provide more precise efficiency corrected cumulants when using the large number of efficiency bins.
	However, it is difficult to know how many efficiency bins are enough. 
	
  * In order to obtain the averaged efficiency for each efficiency bin, we need to consider 
	the variation of the particle yields within the efficiency bin, which means the corrected spectra 
	is necessary. This is especially crucial issue for new data set in 
	new collision energy or collision system where the corrected 
	spectra is not available.
	
  * The calculation cost on the bootstrap increases as ∝ n^M, 
	which indicates that the statistical error estimation will be difficult 
	for many efficiency bins in view of the computing source. 





 §.§ Toy model simulation


Let us discuss more about the first issue related to the approximation 
of efficiency bins above by using a simple toy model.
We generated 50 independent binomial distributions 
for positively and negatively charged particles, 
P^±_i(N^±_i) (i=1,2,...,50). 
The parameters of the binomial distributions are selected to be 
N^+=4, N^-=3 and ε^±=0.8 in Eq. (<ref>). 




Efficiencies for each distribution (ε^±_i) are randomly chosen within   
(0.5,1.0) with uniform distribution.




Particles N^±_i are then randomly sampled with binomial efficiencies ε^±_i 
to define the measured particles n^±_i.
Hereafter, let us suppose the net-particle distribution 
which consists of 50 independent distributions given by
P_ net(N) with N=∑_i=1^50(N^+_i-N^-_i).
We define the m-th order cumulant 
of P_ net(N) as "true" cumulants.
The efficiency correction is performed with 
M efficiency bins by using the corresponding averaged efficiency.  
For instance, let us suppose the efficiency correction with M=10 
efficiency bins. In this case, whole 50 distributions are 
equally divided into 10 sub-bins, 
each sub-bin contains 5 distributions. 
The number of measured particles are counted at each sub-bin 
n^±_ sub,x, (x=1,2,...,10). 
Also the averaged efficiency at each sub-bin is given by 

    ε^±_ sub,x
    	=⟨∑_i=5x^5(x+1)ε^±_iN^±_i
    	/∑_i=5x^5(x+1)N^±_i⟩,

where the bracket represents the event average.
Then n^±_ sub,x and ε^±_ sub,x are substituted 
into Eqs. (<ref>)–(<ref>) 
event-by-event to calculate cumulants
[ Since the notation of electric charge is 
explicitly included in Eq. (<ref>), 
we need following modification for substitution: 
M→2M, a_i=+1  (i≤ M), a_i=-1  (i>M) ].
The efficiency correction has been done with different number of efficiency 
bins for M=1, 2, 5, 10, 25, 50 to check how many efficiency bins are needed 
to obtain the true cumulants.
Figure <ref> shows the corrected 
C_2, C_3 and C_4 as a function 
of the number of efficiency bins, 
where the true value of cumulants are shown in red squares. 
We find that the results with M=50 are consistent with the true cumulants, which is because 50 efficiency bins are assumed in the toy model.
It is also found that using wide efficiency bins is clearly incorrect and 
M=25 efficiency bins is still not enough. 









§ SOLUTION




 §.§ Track-by-track efficiency

























Let us suppose infinite number of efficiency bins M→∞. 
Equation (<ref>) is then written by


    q_(r,s)   =   ∑_i=1^∞ (a_i^r/ε_i^s) n_i
       =   a_1^r/ε_1^sn_1 
    	+ a_2^r/ε_2^sn_2 
    	+ ⋯
    	+ a_i^r/ε_i^sn_i 
    	+ ⋯,

Since we consider that the width of the efficiency bin is now zero, 
each efficiency bin contains up to one particle. 
Thus, the summations for n_i=0 vanish 
(in other words, efficiency bins containing no particles 
don't need to be taken into account), 
and we immediately find that Eq. (<ref>) 
is equivalent to the summation with respect to the 
total number of particles n_ tot in one event, 
which is given by

    q_(r,s) = ∑_j=1^n_ tota_j^r/ε_j^s,

which is connected with Eq. (<ref>) 
via n_ tot = ∑_i=1^M n_i.
It is found that no variable related to the 
efficiency bin appears in Eq. (<ref>). 
What we need to take care of is 
only the track-by-track efficiency, 
which indicates that 
the analytical formula of the efficiency 
with respect to track-wise variables 
can be directly used to determine 
the track-by-track efficiency 
[Dependence of the efficiency 
on event-wise variable 
can be also included if necessary.].
Accordingly, we don't need to estimate 
the averaged efficiency at each efficiency bin, 
so the corrected spectra is no longer 
necessary for the efficiency correction.
The first two difficulties have been solved. 
Hereafter, let us call the correction 
in Eq. (<ref>) "bin-by-bin" method, 
and call the track-wise correction in Eq. (<ref>) 
"track-by-track" method.



 §.§ Method validation

In this sub-section, we employ a toy model 
to demonstrate the validity of the track-by-track efficiency method, 
and also discuss the rest one problem regarding 
how to estimate the statistical errors.
We start from two Gauss distributions, 
one is for positively charged particles, 
and the other is for negatively charged particles.
Figure <ref>–(a) shows the even-by-event correlation 
histogram between positively N^+ and negatively N^- 
charged particles. 
For each particle p_T is allocated according to the 
spectra given by

    f(p_T)=p_Te^-p_T/t  (0<p_T<2),

where t=0.26 and 0.24 for positively and negatively charged 
particles, respectively.
p_T dependent efficiency is assumed to be the convolution 
of the 1st and 2nd polynomial functions 
given by

    ε(p_T) =  a × p_T^2 + b× p_T + c   (0<p_T<1),   
     
    	d× p_T + e   (1<p_T<2),

where (a,b,c,d)=(-0.2,0.4,0.6,0.15,0.65) and (-0.2,0.2,0.65,0.1,0.55) 
for positively and negatively charged particles, respectively.
Each particle is then sampled by the corresponding value of efficiency 
determined in Eq. (<ref>). 
The resulting correlation between measured positively n^+ and negatively 
n^- charged particles is shown in Fig. <ref>–(b), 
and p_T spectra is shown in Fig. <ref>–(d).

Since efficiency changes continuously with respect to p_T, 
it is cumbersome to use the bin-by-bin correction 
formulas in Eq. <ref>. 
Instead, we substitute the value of efficiency for each measured particle
determined by Eq. (<ref>) into Eq. <ref> 
to calculate q_(r,s) at each event. 

The statistical errors can be estimated by the bootstrap method. 
As was mentioned in Sec. <ref>, 
in the case of bin-by-bin correction method 
with M efficiency bins, 
the bootstrap sampling would be performed based on 
the M dimensional histogram.
But now there is no longer the efficiency bin, 
sampling can be simply done based 
on the spectra as follows:

	
  * Resample n_1 and n_2 randomly from Fig. <ref>–(b)
		 
	
  * Allocate p_T for each particle based on the measured spectra in Fig. <ref>–(d)
		 
	
  * Apply the efficiency correction to obtain efficiency corrected cumulants by using the known efficiency curve in Fig. <ref>–(c)
		 
	
  * Repeat <ref>–<ref> with 100 times 
		and take the standard deviation as the statistical error. 
		 

Above procedures are repeated with 100 times independently in order 
to check the validity of the correction and its statistical error. 
Results up to the fourth order are shown in Fig. <ref>. 
It is found that the data points are distributed around the true value, 
so the correction method using track-by-track efficiency works well.
The probability of data points touching the true value 
within the statistical error is shown in the top right of each panel.
We find it comparable with the 1σ nature of the Gaussian, 
which indicates the validity of the bootstrap.
We also checked that the calculation cost only depends 
on the number of particles ∝ n_ tot, 
which is due to the fact that we don't have to 
consider particles which was not measured, 
while the bins containing no particles needed 
to be taken into account in bin-by-bin method.
On the other hand, we can also use the Delta theorem <cit.> to calculate the 
statistical errors of the efficiency corrected cumulants, which will be implemented in the future data analysis. 

It should be noted that one would need to consider how to treat the momentum resolution in real 
experiments. Since we use p_T of individual particles for the efficiency correction, the momentum resolution 
might directly affect the cumulants. This effect can be studied by smearing p_T for each particle with the known value of 
momentum resolution. Finally, we note that bin-by-bin method (even single efficiency bin) 
using the averaged efficiency should also work in this toy model.
This is because only one probability distribution function 
is considered, which indicates we assume that 
the underlying physics is identical 
for whole p_T region for each 
electric charge <cit.>.











§ SUMMARY


In this paper, we pointed out three difficulties 
arising in the experimental application of the 
current efficiency correction method based on the 
arbitrary binning with respect to track-wise variables.
It was shown that those difficulties can be addressed by using the 
track-by-track efficiency in the efficiency correction formula 
based on factorial cumulants. 
We don't need to worry about how many efficiency bins are enough 
to calculate the efficiency corrected cumulants by using the analytical 
parameterization of the efficiency directly.
Thus, the averaged efficiency doesn't need to be estimated and the cumulant analysis can be proceeded without the corrected spectra.
Furthermore, the calculation cost for the statistical error estimation have been significantly reduced. 
Experimentally, single particle efficiency of the detectors 
would depend on transverse momentum, rapidity and azimuthal angle 
in one event. We assume that efficiencies for individual particles 
are independent, which leads to the binomial response function 
of the single particle efficiency.
One thing we have to work hard is to parametrize the single 
particle efficiency as a function of track-wise variables 
with the best precision as long as the computing source allows.
On the other hand, one should also study the possible non-binomial
efficiency effects in the experimental condition. Finally, we emphasize that the method shown in this paper could serve as 
one of the most precise and efficient way for the cumulant efficiency correction with binomial response function.
It will play an important role for the QCD critical point search and can be applied for the cumulant analysis in the future heavy-ion collision experiments, such as 
the BES-II program at RHIC, experiments at FAIR and NICA facilities.



§ ACKNOWLEDGEMENT

This work is supported by the MoST of China 973-Project No. 2015CB856901, 
the National Natural Science Foundation of China under Grants 
(No. 11575069, 11828501, 11890711 and 11861131009), Fundamental Research Funds for the Central Universities No. CCNU19QN054 and China Postdoctoral 
Science Foundation funded project 2018M642878.



