The method we describe here was first developed in [16]. Aside from the theoretical analysis done in [16], we shall also discuss some of the practical considerations when applying this method to the data of Daya Bay.
The key to measuring9Li and8He is the observation that these isotopes are correlated with muons while other events (IBDs originated from the reactors, natural radiations etc.) are not. With this in mind, we try to find the temporal distribution between IBD candidates and the muons preceding them which we shall refer to as time-since-last-muon distribution from now on.
First of all, we consider events that are uncorrelated with muons. Assuming the muons follow Poisson statistics, the inter-arrival time of the muons can be modeled using the exponential distribution with rate parameter Rµ(Tinter-arrival∼ Rµe−Rµt). For uncorrelated signals, using the time reversal property of Poisson process, the time-since-last-muon distribution is simply
funcorr(t) = Rµe−Rµt. (3.2)
Muon induced isotopes, on the other hand, have a different temporal distribution due to their correlations with muons. To derive the distribution for these isotopes, we have to consider cases where the decay windows are intervened by other uncorrelated muons. First of all, let us find the joint distribution of the number of intervening muons and the time-since-last-muon, denoted by fiso(t, n) where n denotes the number of intervening muons as illustrated in figure 3.3. For exponential decay with lifetime τ , Tdcan be modeled by
9Li 13.6
11.8 11.3
7.9
2.8 2.4
9Be
4.7
1.7
8Be+n 3%
1%
2%
16%
30%
49%
β−decay
Figure 3.1: Ground state decays of lithium-9[24]. In Daya Bay, we are mostly concerned with decays into8Be+n for they have similar features to that of IBDs. Note that8Be may further decay into 2 α particles which is not shown in the diagram.
8He 10.7
9.7
5.4
3.2
1.0
8Li
2.5 1%
16%
84%
β−decay
7Li+n 2.0
Figure 3.2: Ground state decays of helium-8.[24] In Daya Bay, we are mostly concerned with decays into7Li+n for they have similar features to that of IBDs.
time
µ0 µ1 µn I
Tn t
Td
...
Figure 3.3: Illustration of time-since-last-muon for muon induced isotopes. In this example, the isotope I is originated from muon µ0 and n uncorrelated muons (µ1,2,...n) intervened this process. Here Td is the decay time of isotope I, Tn is the time interval between µ0 and µnand t is the time-since-last-muon for isotope I. Notice that there can be no other muons between µnand I otherwise µnwould not be the last muon before I.
fTd(t) = 1
τe−τt (3.3)
Since the muons follow Poisson statistics, Tnis the sum of n independent and identically distributed random variables that follow the exponential distribution with rate constant Rµ. In fact, this distribution is widely used in various fields of study and known as the Erlang distribution (see for a proof of this) which has the following form
fTn(t) = Rnµtn−1e−Rµt
(n− 1)! . (3.4)
We are now ready to derive the distribution of t. Notice that t is basically Td− Tnwith an extra condition that no other muons intervened between µnand I (otherwise µnwould not be the last muon before I). This probability distribution can be calculated by
fiso(t, n) = fTd−Tn(t)P(no other muon between µnand I) (3.5)
The marginal distribution fiso(t) is simply
fiso(t) =
Knowing the distributions for both IBDs and muon induced isotopes, the distribution for all IBD candidates can be expressed by
f (t) = αfuncorr(t) + (1− α)fiso(t) (3.13)
where α represents the fraction of uncorrelated events in the IBD candidates. In practice, there may be more than one kind of isotopes that have to be considered (e.g.9Li,8He,12B etc.). In that case a more general form can be used
f (t) = α0funcorr(t) +
isotopesX
i
αifi(t), α0+
isotopesX
i
αi = 1 (3.14)
where fiis the time-since-last-muon distribution for isotope i and αiis the fraction of each isotope in the IBD candidates sample.
Estimation of αican be done using maximum likelihood which has been widely used in the statistics literature for parameter estimation. The various properties and applications in experimental particle physics of maximum likelihood estimation can be found in the statistics chapter of [25]. In our case, given a time-since-last-muon dataset t = {ti}, α ={αi} can be estimated by maximizing the quantity
L(t|α) =Y
i
f (ti|α) (3.15)
which is called the likelihood function. In other words, we want to find ˆα such that
L(t| ˆα) ≥ L(t|α) ∀α ∈ {set of hypotheses} (3.16)
However, care must be taken when applying this model. When muon rates are high, the two distributions can degenerate in the sense that
funcorr(t)≈ fisotope(t) (when Rµ ≫ τ−1)
which would lead to increased uncertainties in our estimations of α. In [16], this effect was estimated to be
σα = 1
√N q
(1 + τ Rµ)2− 1 (3.17)
where N is the total number of events. In fact, this is one of the most crucial problem in
9Li and8He estimation in Daya Bay. Because the muon rates in the near sites were in the order of 10 Hz. A rough estimation using equation 3.17 shows that σα ∼ 0.4%, which
the previously estimated contamination of9Li and8He was about 0.4% in EH1 AD1[17].
In other words, the uncertainty is about 100% were the model applied naively. For this reason, studies were done to include additional physical properties to further constraint the estimations as we shall discuss in the following sections.
3.2.1 An interesting property
Careful readers might notice that during the derivation of the time-since-last-muon distribution, the only assumption we made about the muons was that they follow Poisson statistics which is not related to the actual physical properties of muons. In fact, this assumption can lead to some interesting consequences. First, we can adapt a loose muon selection criteria and include as much muons as possible as long as the non-muon events are independent of muons and follow Poisson statistics. Second, we can divide the muon samples into several mutually exclusive subsets and measure the contaminations of9Li/8He independently and sum them up to get the final result since splitting a Poisson process would result in independent Poisson processes.
In practice, we divide muons by their energies because we know that high energy muons are more likely to produce9Li and8He due to larger cross sections [22] and their rates are low. If we were to simultaneously estimate in both high and low energy muon regions, we would only lose our precision due to the high muon rates from low energy muons. By separating these muons, we can reduce the problem into estimations of the contamination in low energy muon regions.
3.2.2 Practical considerations
In practice, instead of directly maximizing the likelihood function (equation 3.15), we usually solve the problem by minimizing the negative log likelihood
− ln L = −X
i
ln f (ti|α) (3.18)
since the gradient calculations are much easier in the log likelihood formulation which are essential to many numerical optimization algorithms, and both problems share the same solutions due to the monotonicity of the log function. Furthermore, the uncertainty of ˆα
can be estimated, assuming large sample limits, by
( ˆV−1)ij =−∂2lnL
∂αi∂αj (3.19)
where V is the covariance matrix. The corresponding estimated statistical uncertainties for each parameter αiis simplyp
Vˆii.
Notice that the computational resources required to compute eqaution 3.18 scales linearly with the amount of data we have which can be very challenging for experiments such as Daya Bay where millions of events were acquired. A remedy for this problem is to use Poisson binned likelihood estimation where the dataset t is first grouped into several bins where each bin is expected to follow Poisson statistics with expected number of events being
ˆ nb = N
Z
bin
dtf (t|α) (3.20)
where N is the size of the dataset. The log likelihood can then be calculated by
− ln L = −X
b
ln P (nb; ˆnb) (3.21)
=−X
b
(nbln ˆnb− ˆnb− ln Γ(nb+ 1)) (3.22)
where P (n; k) is the Poisson distribution with rate parameter k. The computational complexity now scales with the number of bins which can be controlled depending on the computational resources available.
Statistical and numerical problems of this kind have been studied extensively and various softwares were developed to solve these problems. In the following analysis, we shall useROOT, an object oriented framework for large scale data analysis written in C++ developed at CERN.[26] This framework contains various commonly used statistical and numerical methods in experimental particle physics, and have proven to be effective under such uses.
3.2.3 A naive attempt
Here we give an example on why naively applying the model may lead to unreliable results. The data presented here followed mostly the nGd selection criteria described in the previous chapter with slight modifications which shall be discussed in detail in the next chapter. The muons in this analysis were divided into three regions by the detected number of photoelectrons (p.e.)1, namely the low (p.e. ∈ (3 × 103, 1× 105)), mid (p.e.∈ (1 × 105, 3× 105)) and high (p.e. > 3× 105) regions. Due to the scarcities of cosmogenic isotopes, the fits were done on a site-by-site basis instead of the AD-by-AD measurements. Here we are interested in getting the estimation of the sum of N9Li in all three Eµ regions. However, as one can see from the fits shown in figure 3.4, the uncertainty of N9Lifrom the low Eµregion is so large (2.3× 105) that it alone could ruin the entire estimation.
But let’s not lose our hope, one can see that by dividing the muons into regions, we managed to measure the high Eµ region with much lower uncertainty. The remaining question is how we can improve the resolution of the low/mid Eµ region measurements.
In fact, as we shall discuss in the next section, being able to measure accurately in the high Eµ region had opened up additional possibilities to constrain our low/mid region measurements.