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.

My previous post discussed the mean and variance of photon counts with deadtime. In this post, I describe a model for the energy spectrum that might be measured by a photon counter with perfect pulse height analysis (PHA). Again, my purpose is to gain insight so I will use a highly simplified model. I derive a theoretical formula for the measured spectrum and then a Monte Carlo simulation to validate the model.

I will assume that we use a non-paralyzable photon counter. Recall that this means that the arrival of a photon triggers a deadtime period where subsequent photons are not counted. Photons within the deadtime do not extend it so the counter is enabled immediately after the deadtime for the first photon.

Suppose we use PHA to measure the energy of recorded photons. I will assume that the measured energy is the sum of the energies of the incident photons during the deadtime. This model is widely used as a first approximation for example by Taguchi et al.[1] and Wang et al.[3]. These papers extend the model to use more realistic models of the pulse shape but the simplified model retains many of the fundamental effects of deadtime on the measured spectrum.

If the deadtime is very short and and no additional photons arrive, then only the initial photon energy is measured and the spectrum will be that of the incident photons. With a finite deadtime, the energies of the photons that arrive during the deadtime are added to the initial and the measured spectrum shifts towards higher energies. Since the spectrum is averaged over many photon arrivals, we can express it as a sum over conditional probabilities of 0 or more additional photons in the deadtime. In the following sections, I will first derive the energy spectrum of the sum of *k* photons, then the probability of *k* photons during the deadtime and finally combine them to derive a formula for the measured energy spectrum.

Instead of dealing directly with spectra, the derivation is simplified by using the probability density function (PDF) of the photon energies, which is simply the normalized spectrum
*p*(*E*) = (*S*(*E*))/(*N*_{0}).
where *N*_{0} = ∫*S*(*E*)*dE*. The energies of the individual photons are independent random variables all with a distribution *p*(*E*). To derive the PDF of the sum of two energies, we can use the fact that the PDF of the sum of independent random variables is the convolution of their individual PDF’s. A useful way to derive this is to use the characteristic function, which is the Fourier transform of the PDF,
*C*(*w*) = ⌠⌡*p*(*v*)*e*^{iwv}*dv*.
Inspecting the integral, we can see that it is the expected value ⟨*e*^{iwv}⟩. Now, if the random variables *x* and *y* are independent, then ⟨*xy*⟩ = ⟨*x*⟩⟨*y*⟩. This factoring is also valid for functions of the independent random variables, ⟨*g*(*x*)*h*(*y*)⟩ = ⟨*g*(*x*)⟩⟨*h*(*y*)⟩. If *z* = *x* + *y*, then its characteristic function is
*C*_{z}(*w*) = ⟨*e*^{iw(x + y)}⟩ = ⟨*e*^{iwx}*e*^{iwy}⟩ = ⟨*e*^{iwx}⟩⟨*e*^{iwy}⟩ = *C*_{x}(*w*)*C*_{y}(*w*)
where the last step follows because *x* and *y* are independent. Applying the convolution theorem of Fourier transforms
*p*_{z} = *p*_{x}**p*_{y}
where the symbol * indicates a convolution.

For *k* photons, we apply the theorem multiple times to show that
where *p*_{0}(*E*) = *p*(*E*). To get the energy spectrum, we multiply by *N*_{0} = ∫*S*(*E*)*dE*. Then the spectrum for *k* photons during the deadtime is *S*_{k}(*E*) = *N*_{0}*p*^{(k)}**p*.

I am assuming that the incident counts are a Poisson process. Therefore, the number of counts from the time of the first photon *t*_{0} to *t*_{0} + *τ*, where *τ* is the deadtime is a Poisson random variable with expected value *λ**τ*. The probability of *k* photons during this time is

Combining Eqs. 2↑ and 3↑ with Eq. 1↑,
This formula is implemented with *SpectrumWithPileup*.*m* using the Matlab *conv* function. The number of terms is chosen so the sum of the terms dropped is negligible.

The simulation works in two steps. First, it generates random interarrival times and energies for the incident photons. Next, it processes the arriving photons one-by-one to compute the recorded photon counts and their energies. The recorded energies are computed using the non-paralyzable model as the sum of all the energies of the initial photon and the photons arriving during the deadtime.

The random variables are computed with the inverse transform algorithm. I will use this often so I discuss it now in some detail. I will use the derivation in sec. 5.1 of Ross, **Simulation[2]**. The algorithm is based on the following theorem

LetUbe a uniform(0, 1) random variable. For any continuous cumulative distribution function (CDF)F, the random variable

has distributionF. The notationF^{ − 1}is the inverse function soF^{ − 1}(F(x)) =x.

By definition, the cumulative distribution function of *X* is *a* ≤ *b* as *F*(*a*) ≤ *F*(*b*). Therefore the inequality in the last step of Eq. 6↑ is
*F*_{X}(*x*)
=
*Prob*{*F*(*F*^{ − 1}(*U*)) ≤ *F*(*x*)}
=
*Prob*{*U* ≤ *F*(*x*)}
=
*F*_{U}[*F*(*x*)]
=
*F*(*x*)
where, in the last step, I used the fact that the cumulative distribution function of a uniform random variable is *F*_{U}(*u*) = *u*, 0 ≤ *u* ≤ 1.

(6)
*F*_{X}(*x*)
=
*Prob*{*X* ≤ *x*}
=
*Prob*{*F*^{ − 1}(*U*) < *x*}

where in the last step, I used Eq. 5↑. Since any CDF is monotonically increasing function, we can rewrite any inequality
As an example of the use of this theorem, we can generate random photon interarrival times based on the PDF of the exponential random variable
*p*(*δ**t*) = *λ**e*^{ − λδt}, *δ**t* ≥ 0
The cumulative distribution function is *F*(*δ**t*) = *λ*∫^{δt}_{0}*e*^{ − λt}*dt* = 1 − *e*^{ − λδt} so its inverse is
*F*^{ − 1}(*u*) = − (1)/(*λ*)log(1 − *u*).
Since 1 − *u* is also a uniform(0, 1) random variable, we can generate the interarrival times using the Matlab *rand* function to generate *ntimes* values
*dt*_*pulse* = ( − 1 ⁄ *lambda*)**log*(*rand*(*ntimes*, 1)).

With the energy spectrum, we do not have an analytical inverse function but we can use the Matlab *cumsum* and *interp*1 functions. See the code for *SpectrumCumulativeInverse*.*m*.

The *EgysWithPileup*.*m* function uses these formulas to generate random photon interarrival times and the corresponding random photon energies. Then it generates the recorded energies using the following Matlab function. Note the similarity of the function to the *CountPulses* function of my previous post.

The code for this post has a function to carry out the Monte Carlo simulation. The code computes the random recorded energies with deadtime. Fig. 1↓ shows a comparison of a histogram of the Monte Carlo results and the spectrum computed using Eq. 1↑ for two values of the mean incident counts per second times deadtime *η* = *λ**τ* parameter. The histogram bins are computed in my *HistOptBin*.*m* function using the Freedman-Diaconis rule, which is based on number of data points and the inter-quartile range.

Deadtime has a strong effect on the measured spectrum smearing it to higher energies than are present in the input spectrum. The data in Fig. 1↑ were computed with parameters *η* = *λ**τ* = 0.6 and 0.1 This is the mean number of incident photons during the deadtime. With high count rates, it is difficult to achieve a low value. For example, if the incident rate is 10^{8} counts per second, then an *η* ≤ 0.6 implies the deadtime is less than 6 nanoseconds. Taguchi et al.[1] report measured deadtimes of approximately 100 nanoseconds with year 2011 state of the art detectors.

The situation is not as bleak as these numbers indicate because the 10^{8} incident rate is only in the “air” data outside the patient. With the collimation in CT systems, the body transmission in the center can be less than 0.01. In this case, the mean counts per deadtime woud be *η* = 10^{6}*100*10^{ − 9} = 0.1. As shown in the right panel of Fig. 1↑, the distortion is much smaller in this case. Nevertheless, higher rates are present at the edges of the patient and we need to develop algorithms that can handle the distorted data in these regions.

Last edited Oct. 08, 2011

© 2011 by Aprend Technology and Robert E. Alvarez

Linking is allowed but reposting or mirroring is expressly forbidden.

[1] : “Modeling the performance of a photon counting x-ray detector for CT: Energy response and pulse pileup effects”, *Med. Phys.*, pp. 1089—1102, 2011.

[2] : *Simulation, Fourth Edition*. Academic Press, 2006.

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