If you have trouble viewing this, try the pdf of this post. You can download the code used to produce the figures in this post.

In previous posts, I discussed the mean and variance and the energy spectrum of photon counting with deadtime. In this post, I will describe the statistics of pulse height analysis (PHA) data as a function of the deadtime of the detector. I will analyze the idealized case with perfect energy bins with zero transition width and no overlap and no added electronic noise. With these assumptions and no deadtime, the number of counts in each bin is Poisson distributed with a mean value equal to the number of incident photons and the data in different bins are independent. With deadtime, the PHA data mean and variance are smaller than those with no deadtime. In addition, the data in different bins become negatively correlated.

In my next post, I will describe a Monte Carlo simulation to validate the formulas derived here.

The derivation here is taken from a recent paper by Wang et al.[2]. The approach is first to derive statistics with a fixed number of recorded counts and then to use conditional expectation, variance, and covariance to derive the statistics with random numbers of recorded counts. By fixing the recorded counts, the counts for the PHA bins become multinomial random vectors. So, first I will derive the basic properties a multinomial random vector. Then, I will summarize the formulas for conditional statistics and apply them to the PHA count data.

The multinomial distribution is a generalization of the binomial random variable that is discussed in most probability textbooks. You will recall that a binomial is the number of successes in *n* independent trials if the probability of success in each trial is *p*. With a multinomial, there are more than two possible outcomes of each trial ( see Ross [1] Section 6.1). If there are *n* trials with *r* possible outcomes with probabilities *p*_{r} *r* = 1…*r* and *X*_{j}, *j* = 1…*r* is the number of successes for outcome *j* , then the probability mass function of the joint distribution of the *X*_{j} is

(1) *Prob*(*X*_{1} = *k*_{1}, *X*_{2} = *k*_{2, }, …, *X*_{r} = *k*_{r}) = (*n*!)/(*k*_{1}!*k*_{2}!⋯*k*_{r}!)*p*^{k1}_{1}*p*^{k2}_{2}⋯*p*^{kr}_{r}

I will derive the moments of the multinomial distribution using the moment generating function, which, for a joint distribution, is defined to be (see Ross [1] Section 7.6)
*M*_{X}(*t*_{1}, …, *t*_{r}) = ⟨*e*^{X1t1 + ⋯ + Xrtr}⟩

In my notation ⟨⋅⟩ is the expected value. This function can be derived by using the definition of expected value with the probability mass function in Eq. 1↑

(2)
⟨*e*^{X1t1 + ⋯ + Xrtr}⟩
=
^{ }⎲⎳_{{ki:⎲⎳ki = n}}(*n*!)/(*k*_{1}!*k*_{2}!⋯*k*_{r}!)*p*^{k1}_{1}*p*^{k2}_{2}⋯*p*^{kr}_{r}*e*^{k1t1 + ⋯ + krtr}
=
^{ }⎲⎳_{{ki:⎲⎳ki = n}}(*n*!)/(*k*_{1}!*k*_{2}!⋯*k*_{r}!)(*p*_{1}*e*^{t1})^{k1}(*p*_{2}*e*^{t2})^{k2}⋯(*p*_{r}*e*^{tr})^{kr}

The multinomial theorem from elementary algebra gives the expansion of a multinomial to the n’th power

(3) (*x*_{1} + *x*_{2} + ⋯ + *x*_{r})^{n} = ^{ }⎲⎳_{{ki:⎲⎳ki = n}}(*n*!)/(*k*_{1}!*k*_{2}!⋯*k*_{r}!)*x*^{k1}_{1}*x*^{k2}_{2}⋯*x*^{kr}_{r}

Comparing this to Eq. 2↑, we see that by substituting *x*_{k} = *p*_{k}*e*^{tk}

(4) *M*_{X}(*t*_{1}, …, *t*_{r}) = ⟨*e*^{X1t1 + ⋯ + Xrtr}⟩ = (*p*_{1}*e*^{t1} + *p*_{2}*e*^{t2} + ⋯ + *p*_{r}*e*^{tr})^{n}

We can use the moment generating function in Eq. 4↑ to derive the moments using the general relation

⟨*X*^{j1}_{1}⋯*X*^{jr}_{r}⟩ = ⎡⎣(∂*M*_{X}(*t*_{1}, …, *t*_{r}))/(∂*t*^{j1}_{1}⋯∂*t*^{jr}_{r})⎤⎦_{t1, …, tr = 0}

Taking the first derivative, the expected value of any component of the multinomial is

⟨*X*_{k}⟩ = [*n*(*p*_{1}*e*^{t1} + *p*_{2}*e*^{t2} + ⋯ + *p*_{r}*e*^{tr})^{n − 1}*p*_{k}*e*^{tk}]_{t1, …, tr = 0} = *np*_{k}

since
^{r}⎲⎳_{k = 1}*p*_{k} = 1.

Taking another derivative, the second moment is ⟨*X*^{2}_{k}⟩ = *n*(*n* − 1)*p*^{2}_{k} + *np*_{k} so the variance is
*Var*(*X*_{k}) = ⟨*X*^{2}_{k}⟩ − ⟨*X*_{k}⟩^{2} = *np*_{k}(1 − *p*_{k})

Notice that the individual counts have the same mean and variance as a binomial random variable. This makes sense since in any trial, a particular case will occur with probability *p*_{k}.

Taking derivatives with respect to *t*_{j} and *t*_{k}, the covariance for *j* ≠ *k* is
*Cov*(*X*_{j}, *X*_{k}) = − *np*_{j}*p*_{k}

Of course, if *j* = *k*, *Cov*(*X*_{k}, *X*_{k}) = *Var*(*X*_{k}).

Ross[1] discusses conditional expectation in Section 7.4. He proves the following result

⟨*X*⟩ = ⟨⟨*X*|*Y*⟩⟩_{Y}

In addition, I will give a proof for the conditional covariance formula (see Wikipedia)

This can be proved by starting with the computational formula for covariance
*cov*(*X*, *Y*) = ⟨*XY*⟩ − ⟨*X*⟩⟨*Y*⟩

Rewriting the right hand side by using conditional expectation with the variable *Z*,

Using the covariance computational formula again, ⟨*XY*|*Z*⟩ = *cov*(*X*, *Y*|*Z*) + ⟨*X*|*Z*⟩⟨*Y*|*Z*⟩. Substituting in the first term on the right hand side of (6↑)

(7)
*cov*(*X*, *Y*)
=
⟨*cov*(*X*, *Y*|*Z*) + ⟨*X*|*Z*⟩⟨*Y*|*Z*⟩⟩ − ⟨*X*|*Z*⟩⟨*Y*|*Z*⟩
=
⟨*cov*(*X*, *Y*|*Z*)⟩ + ⟨⟨*X*|*Z*⟩⟨*Y*|*Z*⟩⟩ − ⟨*X*|*Z*⟩⟨*Y*|*Z*⟩

where I have used the fact that the expectation of a sum is the sum of the expectations. The last two terms of (7↑) are the computational formula for *cov*(⟨*X*|*Z*⟩, ⟨*Y*|*Z*⟩), so we have derived the conditional covariance formula, Eq. 5↑.

Since *var*(*X*) = *cov*(*X*, *X*), we can use (5↑) to show that the formula for conditional variance is

Now we are ready to derive the statistics of PHA data with deadtime. The derivation here will follow that in Wang et al.[2]. As I discussed in my blog post, the spectrum of the measured energies with deadtime is

where *s*(*E*) = (*S*(*E*))/(*N*_{0}), *S*(*E*) is the incident energy spectrum, and *N*_{0} = ∫*S*(*E*)*dE*. Each term in the sum corresponds to *k* additional photons arriving during the dead time and, by definition, *s*^{(0)}**s* = *s*.

With PHA, we group the measured energies into *r* bins, corresponding to energies from [*E*_{k − 1}:*E*_{k}, *k* = 1…*r*]. Suppose we have a fixed number of recorded counts, *M*. Then, since the photons have random energies, each bin *M*_{k}|*M*, *k* = 1…*r* counts the number of photons with energies in the bin energy range. The probability for each bin is
*p*_{k} = (^{Ek}⌠⌡_{Ek − 1}*S*_{deadtime}(*E*)*dE*)/(⌠⌡*S*_{deadtime}(*E*)*dE*)

and the counts are components of a multinomial random vector. By the results above, the expected value and the variance of the counts in a bin for fixed *M* are
*var*(*M*_{k}|*M*) = *Mp*_{k}(1 − *p*_{k}).

⟨*M*_{k}|*M*⟩ = *Mp*_{k}

Applying the conditional expectation formula, the mean value for each bin is

(10)
⟨*M*_{k}⟩
=
⟨⟨*M*_{k}|*M*⟩⟩
=
⟨*Mp*_{k}⟩
=
⟨*M*⟩*p*_{k}

We can apply the laws of conditional variance Eq. 8↑ and covariance Eq. 5↑ to derive the statistics of the PHA data. First the variance

(11)
*var*(*M*_{k})
=
⟨*var*(*M*_{k}|*M*)⟩ + *var*(⟨*M*_{k}|*M*⟩)
=
⟨*Mp*_{k}(1 − *p*_{k})⟩ + *var*(*Mp*_{k})
=
*p*_{k}(1 − *p*_{k})⟨*M*⟩ + *p*^{2}_{k}*var*(*M*)
=
⟨*M*⟩*p*_{k} + (*var*(*M*) − ⟨*M*⟩)*p*^{2}_{k}

The covariance can be derived similarly

(12)
*cov*(*M*_{j}, *M*_{k})
=
⟨*cov*(*M*_{j}, *M*_{k}|*M*)⟩ + *cov*(⟨*M*_{j}|*M*⟩, ⟨*M*_{k}|*M*⟩)
=
− *p*_{j}*p*_{k}⟨*M*⟩ + *p*_{j}*p*_{k}*var*(*M*)
=
*p*_{j}*p*_{k}(*var*(*M*) − ⟨*M*⟩)

The expected value in Eq. 10↑ is straightforward but *p*_{k} is the fraction of the distorted spectrum with non-zero deadtime, Eq. 9↑, in each energy bin. In a typical x-ray imaging system, this spectrum and therefore the fractions can change markedly since the incident count rate can vary by a factor of 100 from air to the interior of the object.

If the deadtime is zero, then the total recorded counts *M* are Poisson distributed with mean ⟨*M*⟩ = *λ**T* that is equal to the variance *var*(*M*) = *λ**T*. Therefore, by Eq. 11↑, the variance of the counts in each bin is also equal to the expected value. With non-zero deadtime, my previous post shows that the expected value and variance for large counts are
*var*(*M*) = (*λ**T*)/((1 + *λ**τ*)^{3})

⟨*M*⟩ = (*λ**T*)/(1 + *λ**τ*)

where *λ* is the incident count rate, T is the integration time, and *τ* is the deadtime. Note that in the formula for the spectrum with deadtime, Eq. 9↑, above, *N*_{0} = *λ**T*. In this case, *var*(*M*) < ⟨*M*⟩, so from Eq. 11↑ the variance of the PHA bin counts is also less than expected from Poisson statistics, *var*(*M*_{k}) < ⟨*M*⟩*p*_{k}. This makes sense since the photons that arrive during the deadtimes do not change the counts so the variance is reduced.

Also, Eq. 12↑ shows that with zero deadtime so that *var*(*M*) = ⟨*M*⟩, the covariance is 0. With non-zero deadtime, *var*(*M*) < ⟨*M*⟩ and the covariance is negative.

Last edited Oct. 17, 2011

© 2011 by Aprend Technology and Robert E. Alvarez

Linking is allowed but reposting or mirroring is expressly forbidden.

[1] : *A First Course in Probability*. Prentice Hall College Div, 1997.

[2] : “Pulse pileup statistics for energy discriminating photon counting x-ray detectors”, *Medical Physics*, pp. 4265—4275, 2011.