User:Wwheaton/test10

Gamma-ray detector matrices and invertibility -- DRAFT

Wm. A. Wheaton From LSU, 07:27 1 August 2010 (UTC)

Factoring the forward spectrometer matrix G = R*C*H into three factors, the efficiency matrix H (diagonal), Compton matrix C (upper triangular), and resolution matrix R (diagonally dominant, with elements along diagonal roughly like integrals of Gaussian resolution function over energy channel width. Matrix multiplication is "*" here, with energy input on right, pulse height output to the left.  This matrix turns a photon spectrum binned in (possibly narrow) energy channels, first factor into an interaction spectrum, then into an energy loss spectrum in the C matrix, and finally into a PHA spectrum in the output resolution matrix R. Each of these can be factored in turn into SVD matrices   U*S*VT, but the first two (H & C, proceeding from right to left) will be non-singular and well-conditioned, as long as the channel widths are constant, from input to output, so both H & C are square. The matrix R will generally be non-singular if the number of input energy channels is ≤ the the number of output PHA channels, and will be fairly well-conditioned if the the Gaussian sigma is less than the PHA channel width.

In this case, using the SVD theorem on R = UR SR VRT, we may invert G to obtain:

G-1 = {H-1 C-1 VR} SR URT,

where H and C,  are guaranteed to be invertible and well-conditioned, VR is orthogonal, and URT has orthogonal rows. leaving SR as the only possible source of ill-conditioning in the problem. Thus we focus our attention on SR, which is diagonal, but with possibly some 0 or very small diagonal elements. Note in particular that the determinant of SR will be the product of its diagonal elements.

The idea here is to match the output channel widths to the detector resolution matrix R so that the total response matrix G is invertible.

So it seems to me that the strategy is to adjust all the photon energy channel widths so that the columns of the R matrix are all essentially the same, ie, 0,0,0,...0,eps,d,eps,0,0,....,0.where eps+d+eps = 1.000 pretty close, all the d elements are on the diagonal (& equal) and there is a little hanky-panky at the uper and lower edges. In other words, keep the channel a fixed number of sigma wide, so the determinant will come out something like sigma^J, where R is JXJ.

And what I really want to do here is understand the SVD behavior of a nearly diagonal matrix with diagonal elements like d(x) = erf(x)-erf(-x) = 2*erf(x), or something like that....

Of course this means the input channels have to be adjusted to the detector resolution, and the output PHA channels have to be adjusted somehow to fit. I could use very narrow PHA channels, and then sum over them to get output channels any size I want. But how should I do it with data from a real PHA, after the fact, from a real instrument? Evidently the answer is to adjust the energy width of the columns, as those I can set any way I want.

Note that the nnls non-negativity constraint has not been applied yet, but will be, and will make the results still more robust, especially if there is a background subtraction (which, if done carelessly, can give negative fluxes).