About Poisson noise#

In order for Poisson noise to be correctly added to the variance map, a few important details must be understood first.

Theoretical considerations#

We assume that for each pixel \((x, y)\) on the detector, there is a signal \(S(x, y)\) which follows a Poisson distribution with an average value (and variance) \(\lambda (x,y)\). We also assume that the signal is high enough so that \(S\) is a good approximation of \(\lambda\). This implies that \(Var [ S ] \approx S\).

In practice, the variance is not just due to Poisson noise but also to other sources of noise (e.g. readout) which we combine into an additional term \(\sigma^2\). The total variance is therefore:

\[\sigma^2_{\rm tot} (x, y) = S (x, y) + \sigma^2,\]

where \(S\) and \(\sigma^2\) must be in units of counts or \(e^-\).

In the expression above, \(S\) represents the signal measured in a pixel without convolution. With convolution, the signal becomes

\[S' (x, y) = \sum_{x', y'} {\rm PSF} (x' - x, y' - y) S(x', y').\]

It can be shown (e.g. see Appendix D of Mercier W. PhDT) that in such a case and when neglecting covariance between pixels, the variance of \(S'\) writes

\[{\rm Var} [S'] (x, y) = \sum_{x', y'} {\rm PSF}^2 (x' - x, y' - y) \sigma^2_{\rm tot} (x', y'),\]

which can be rewritten as

\[{\rm Var} [S'] (x, y) = \tilde \sigma^2 (x, y) + \sum_{x', y'} {\rm PSF}^2 (x' - x, y' - y) {\rm Var} [S] (x', y'),\]

where \(\tilde \sigma^2 (x, y)\) is the variance map without Poisson noise (but convolved, hereafter background variance). Finally, remembering that \(Var [ S ] \approx S\), we can simplify the expression to

\[{\rm Var} [S'] (x, y) = \tilde \sigma^2 (x, y) + \tilde S (x, y),\]

where \(\tilde S\) is the input image convolved by the square of the PSF (and still in units of counts or \(e^-\)).

Practical implementation#

If Poisson noise is missing from the input variance map, this code will assume that the variance map provided by the user corresponds to \(\tilde \sigma^2\). To add Poisson noise, the first step is to provide the name of a file that contains \(\tilde S\), that is the data convolved by the square of the PSF, when creating the Filter objects:

filt_435w = SED.Filter('F435W', 'data_F435W.fits', 'var_F435W.fits', 25.68,
                       file2 = 'data_conv2_F435W.fits'
                      )

Dealing with units#

Data in \(e^-\,s^{-1}\)#

Because this code was developed with HST images in mind, it assumes that the input data and the data convolved by the square of the PSF are in \(e^-\,s^{-1}\), and that the variance maps are in \(e^-\,s^{-2}\). The code will read from the input image the FITS header keyword 'TEXPTIME' which will be used to convert the variance to the right unit. The conversion is done with the following reasoning:

  • If \(\tilde \sigma^2\) is the background variance in \(e^-\,s^{-2}\), then its value in \(e^-\) must be given by \(\tilde \sigma^2 \times T_e^2\), where \(T_e\) is the exposure time.

  • If \(\tilde S\) is the signal in \(e^-\,s^{-1}\), then its value in \(e^-\) must be given by \(\tilde S \times T_e\).

  • If \({\rm Var} [S']\) is the total variance in \(e^-\,s^{-2}\), then its value in \(e^-\) must be given by \({\rm Var} [S'] \times T_e^2\).

Therefore, using the expression for \({\rm Var} [S']\) and the conversions between \(e^-\) and \(e^-\,s^{-2}\) given above, we find that the total variance in \(e^-\,s^{-2}\) writes

\[{\rm Var} [S'] (x, y) = \tilde \sigma^2 (x, y) + \tilde S (x, y) / T_e,\]

where \(\tilde \sigma^2\) is the background variance map in \(e^-\,s^{-2}\) and \(\tilde S\) is the data in \(e^-\,s^{-1}\).

Data in other units#

The methods setCode() or genTable() can be called with an extra argument to divide the value of the exposure time. For instance

flist.genTable(texpFac=10)

will divide the exposure time by a factor of 10 when computing the Poisson noise contribution to the total variance. If texpFac is equal to 0, no Poisson noise will be added.

This argument can be used to properly handle Poisson noise with data that are not in \(e^-\,s^{-1}\). Indeed, let us assume that the data are in an arbitrary unit \(u\), which can be converted to \(e^-\,s^{-1}\) with a scale factor \(f\) (i.e. \(S [u] = f \times S [e^-\,s^{-1}]\)). Then,

  • If \(\tilde \sigma^2\) is the background variance in unit \(u^2\), its value in \(e^-\) must be given by \(\tilde \sigma^2 \times T_e^2 / f^2\)

  • If \(\tilde S\) is the signal in unit \(u\), then its value in \(e^-\) must be given by \(\tilde S \times T_e / f\).

  • If \({\rm Var} [S']\) is the total variance in unit \(u^2\), then its value in \(e^-\) must be given by \({\rm Var} [S'] \times T_e^2 / f^2\).

Following the same logic as in the previous section, the total variance in unit \(u^2\) writes

\[{\rm Var} [S'] (x, y) = \tilde \sigma^2 (x, y) + f \times \tilde S (x, y) / T_e,\]

where \(\tilde \sigma^2\) is the background variance map in unit \(u^2\) and \(\tilde S\) is the input data convolved by the square of the PSF in unit \(u\).

Thus, the keyword argument texpFac can be used as a conversion factor to properly compute the contribution of the Poisson noise to the total variance.

Data in \(e^-\)#

If the data are in electrons, then no unit conversion is required. The simplest solution is therefore to read the exposure time from the FITS file header and provide it as a scale factor:

from astropy.io import fits

with fits.open('data_F435W.fits') as hdul:
    texp = hdul[0].header['TEXPTIME']

flist.genTable(texpFac=texp)