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 “SNR with pileup ...” paper[1] presented a set of theoretical formulas for the noise of NQ and PHA detectors in Tables I and II. I have discussed the individual formulas and Monte Carlo simulations of their validity in past posts in this series. In Section 2.K of the paper, I presented an overall test of the formulas that compared the A-vector component variances with a Monte Carlo simulation of the random detector data processed with a maximum likelihood estimator (MLE). In this post, I expand the discussion in the paper and present code to reproduce Fig. 2.

Fig. 2 plots the A-vector variances as a function of the number of photons incident on the detector for three values of *η*. The theoretical A-vector components noise variances were taken from the diagonal elements of the Cramèr-Rao lower bound (CRLB) computed with the functions CRLB_NQ_with_pileup.m and CRLB_PHA_with_pileup.m. These functions are “cleaned up” versions of the functions used in previous posts in this series. They use the first difference approach to computing the CRLB with pileup described in this post.

The theoretical variances were compared with the sample variances of a set of A-vector values computed from random samples of detector data. The code to generate the random samples of NQ and PHA data is shown in the text box. It is the same as used in a previous post.

The next step is to calculate A-vectors from the random data with a linear MLE. The linear MLE, discussed in this post, is
**A**̂_{MLE} = [(**M**^{T}C^{ − 1}_{L}M)^{ − 1}** M**^{T}C^{ − 1}_{L}]**L**

Comparing this with the constant covariance approximation to the CRLB, discussed in this post,
**CRLB**_{const − cov} = [**M****C**^{ − 1}_{L}M]^{ − 1},

we can see that the linear estimator is
**A**̂_{MLE} = [CRLB_{const − cov}**M**^{T}C^{ − 1}_{L}]L.

This is the formula implemented in the code.

The MLE were then applied to the logarithm of the random detector data as as shown in the following text boxes for the NQ detector. The **L** data are the negative of the logarithm of the measurements divided by the measurements with no object. For the NQ detector, these “air” values are as as shown. Note that the expected value of Q is ⟨*Q*⟩ = ⟨*N*⟩⟨*E*⟩

Finally, the bootstrap algorithm was used to compute the mean value and the standard deviation of the mean. The standard deviation can be used to display error bars on the Monte Carlo data, although this was not done in the example since they were too small to display.

The PHA results were computed with the same approach.

The Monte Carlo and theoretical data were plotted by the code in the post code file. The results are shown in Fig. 1↓ for the variances vs. the number of photons and in Fig. 2↓ for the variance at fixed number of photons vs *η*.

The Monte Carlo values in Fig. 1↑ and 2↑ are close to the theoretical formulas. These results validate but, of course, do not mathematically prove the formulas. Notice that the NQ variance increases fractionally more than the four bin PHA as *η* increases.

Last edited Feb 10, 2015

Copyright © 2015 by Robert E. Alvarez

Linking is allowed but reposting or mirroring is expressly forbidden.

[1] : “Signal to noise ratio of energy selective x-ray photon counting systems with pileup”, *Med. Phys.*, pp. -, 2014.