aprendtech.com >> blog >> this post

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.

Deadtime-3 Pulse Height Analysis Theory

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

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₁ = k₁, X₂ = k_2,, …, X_r = k_r) = (n!)/(k₁!k₂!⋯k_r!)p^k1₁p^k₂₂⋯p^k_r_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₁, …, t_r) = ⟨e^{X₁t₁ + ⋯ + X_rt_r}⟩

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^{X₁t₁ + ⋯ + X_rt_r}⟩ = ⎲⎳_{{k_i:⎲⎳k_i = n}}(n!)/(k₁!k₂!⋯k_r!)p^k1₁p^k₂₂⋯p^k_r_re^{k₁t₁ + ⋯ + k_rt_r} = ⎲⎳_{{k_i:⎲⎳k_i = n}}(n!)/(k₁!k₂!⋯k_r!)(p₁e^t1)^k₁(p₂e^t2)^k₂⋯(p_re^t_r)^k_r

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

(3) (x₁ + x₂ + ⋯ + x_r)ⁿ = ⎲⎳_{{k_i:⎲⎳k_i = n}}(n!)/(k₁!k₂!⋯k_r!)x^k1₁x^k₂₂⋯x^k_r_r

Comparing this to Eq. 2↑, we see that by substituting x_k = p_ke^t_k

(4) M_X(t₁, …, t_r) = ⟨e^{X₁t₁ + ⋯ + X_rt_r}⟩ = (p₁e^t₁ + p₂e^t₂ + ⋯ + p_re^t_r)ⁿ

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

⟨X^j₁₁⋯X^j_r_r⟩ = ⎡⎣(∂M_X(t₁, …, t_r))/(∂t^j₁₁⋯∂t^j_r_r)⎤⎦_{t₁, …, t_r = 0}

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

⟨X_k⟩ = [n(p₁e^t₁ + p₂e^t₂ + ⋯ + p_re^t_r)^n − 1p_ke^t_k]_{t₁, …, t_r = 0} = np_k

since

^r⎲⎳_k = 1p_k = 1.

Taking another derivative, the second moment is ⟨X²_k⟩ = n(n − 1)p²_k + np_k so the variance is

Var(X_k) = ⟨X²_k⟩ − ⟨X_k⟩² = 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_jp_k

Of course, if j = k, Cov(X_k, X_k) = Var(X_k).

Conditional expectation and variance

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)

(5) cov(X, Y) = ⟨cov(X, Y|Z)⟩ + cov(⟨X|Z⟩, ⟨Y|Z⟩)

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,

(6) cov(X, Y) = ⟨⟨XY|Z⟩⟩ − ⟨X|Z⟩⟨Y|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

(8) var(X) = ⟨var(X|Y)⟩ + var(⟨X|Y⟩)

Statistics of PHA data with deadtime

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

(9) S_deadtime(E) = N₀^∞⎲⎳_k = 0((λτ)^k)/(k!)e^− λτ(s^(k)*s)

where s(E) = (S(E))/(N₀), S(E) is the incident energy spectrum, and N₀ = ∫S(E)dE. Each term in the sum corresponds to k additional photons arriving during the dead time and, by definition, s⁽⁰⁾*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 = (^E_k⌠⌡_{E_k − 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

⟨M_k|M⟩ = Mp_k

var(M_k|M) = Mp_k(1 − p_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²_kvar(M) = ⟨M⟩p_k + (var(M) − ⟨M⟩)p²_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_jp_k⟨M⟩ + p_jp_kvar(M) = p_jp_k(var(M) − ⟨M⟩)

Discussion

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

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

var(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₀ = λ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

Linking is allowed but reposting or mirroring is expressly forbidden.

References

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

[2] Adam S Wang, Daniel Harrison, Vladimir Lobastov, J Eric Tkaczyk: “Pulse pileup statistics for energy discriminating photon counting x-ray detectors”, Medical Physics, pp. 4265—4275, 2011.