diff --git a/harnessed_jobs/read_noise_raft/v0/producer_read_noise_raft.py b/harnessed_jobs/read_noise_raft/v0/producer_read_noise_raft.py index 83d385b..c091491 100755 --- a/harnessed_jobs/read_noise_raft/v0/producer_read_noise_raft.py +++ b/harnessed_jobs/read_noise_raft/v0/producer_read_noise_raft.py @@ -3,11 +3,14 @@ Producer script for raft-level read noise analysis. """ from __future__ import print_function +import matplotlib.pyplot as plt import lsst.eotest.sensor as sensorTest import siteUtils import eotestUtils +from correlated_noise import correlated_noise from multiprocessor_execution import sensor_analyses + def run_read_noise_task(sensor_id): file_prefix = '%s_%s' % (sensor_id, siteUtils.getRunNumber()) bias_files = siteUtils.dependency_glob('S*/%s_fe55_fe55_*.fits' % sensor_id, @@ -26,5 +29,11 @@ def run_read_noise_task(sensor_id): task.run(sensor_id, bias_files, gains, system_noise=system_noise, mask_files=mask_files, use_overscan=True) + # Compute amp-amp correlated noise. + bias_stats, corr_fig, _ = correlated_noise(bias_files, target=0, + make_plots=True, title=sensor_id) + plt.figure(corr_fig.number) + plt.savefig('%s_correlated_noise.png' % file_prefix) + if __name__ == '__main__': sensor_analyses(run_read_noise_task) diff --git a/python/correlated_noise.py b/python/correlated_noise.py new file mode 100644 index 0000000..1f81a1c --- /dev/null +++ b/python/correlated_noise.py @@ -0,0 +1,177 @@ +import glob +from collections import namedtuple, defaultdict +import numpy as np +import scipy +import astropy.io.fits as fits +import matplotlib.pyplot as plt +import lsst.eotest.image_utils as imutils +import lsst.eotest.sensor as sensorTest + +plt.rcParams['xtick.labelsize'] = 'x-small' +plt.rcParams['ytick.labelsize'] = 'x-small' + + +def get_oscan_indices(target_file): + "Return the pixel indices of the overscan region." + amp_geom = sensorTest.makeAmplifierGeometry(target_file) + bbox = amp_geom.serial_overscan + return bbox.getMinY(), bbox.getMaxY(), bbox.getMinX(), bbox.getMaxX() + + +def get_overscans(infile, oscan_indices=None): + "Return the overscan data as numpy arrays." + if oscan_indices is None: + y0, y1, x0, x1 = get_oscan_indices(infile) + else: + y0, y1, x0, x1 = oscan_indices + ccd = fits.open(infile) + overscans = dict() + for amp in imutils.allAmps(): + overscans[amp] = ccd[amp].data[y0:y1, x0:x1] + return overscans + + +def get_mean_overscans(infiles, oscan_indices=None): + "Compute the mean of the overscans of the list of files" + if oscan_indices is None: + y0, y1, x0, x1 = get_oscan_indices(infiles[0]) + else: + y0, y1, x0, x1 = oscan_indices + mean_overscans = defaultdict(list) + for infile in infiles: + ccd = fits.open(infile) + for amp in imutils.allAmps(): + mean_overscans[amp].append(ccd[amp].data[y0:y1, x0:x1]) + for amp, images in mean_overscans.items(): + mean_overscans[amp] = sum(images)/float(len(images)) + return mean_overscans + + +BiasStats = namedtuple('BiasStats', + 'noise_orig noise_corr corr_factor bias_oscan'.split()) + + +def correlated_noise(bias_files, target=0, make_plots=False, plot_corr=True, + figsize=(8, 8), title=''): + """ + Compute the correlated noise statistics for the overscan regions + of the list of files, optionally making plots of the distributions. + + Parameters + ---------- + bias_files: list + List of bias files to analyze. This list must have at least as many + files as the target file index + 1. + target: int + Bias frame to compare to the mean biases constructed from the + remaining files. + make_plots: bool [False] + Flag to determine if the png plots will be generated. + plot_corr: bool [True] + Flag to plot the histograms of correlation-corrected pixel + values. If False, then plot histograms of the uncorrected pixel + values. + figsize: tuple [(8, 8)] + Figure size (in inches) of 4x4 grid of correlation plots. + title: str [''] + Title of 4x4 grid. + + Returns + ------- + (dict, figure, figure): tuple of results and matplotlib + figures. The first item is a dict of BiasStats objects, + + BiasStats = namedtuple('BiasStats', \ + 'noise_orig noise_corr corr_factor bias_oscan'.split()) + + that contain the results for each amplifier. + """ + f1, f2 = None, None + if make_plots: + f1, ax1 = plt.subplots(4, 4, figsize=figsize) + ax1 = {amp: subplot for amp, subplot in zip(imutils.allAmps(), + ax1.flatten())} + f2, ax2 = plt.subplots(4, 4, figsize=figsize) + ax2 = {amp: subplot for amp, subplot in zip(imutils.allAmps(), + ax2.flatten())} + + # Extract the target filename and omit it from the list of bias files. + target_file = bias_files.pop(target) + + # Get the target frame overscans. + bias_oscans = get_overscans(target_file) + oscan_shape = bias_oscans[1].shape + + # Construct the mean bias overscans from the remaining files. + mean_oscans = get_mean_overscans(bias_files) + + # Compute the mean values of the mean bias overscans. + mean_oscan_values \ + = {amp: np.mean(oscan) for amp, oscan in mean_oscans.items()} + + # Loop over amps in target frame and compute statistics. + bias_stats = dict() + for amp in bias_oscans: + # Loop over other amps and construct the mean image of the + # bias-subtracted overscans. Require included amps to have + # (unsubtracted) overscans with 4 < stdev < 25 rms ADU. + reduced_mean_oscan = np.zeros(oscan_shape) + num_oscan = 0 + for oamp, oscan in bias_oscans.items(): + if oamp == amp or not (4. < np.std(oscan) < 25): + continue + reduced_mean_oscan += (oscan - mean_oscans[oamp]) + num_oscan += 1 + reduced_mean_oscan -= np.mean(reduced_mean_oscan) + reduced_mean_oscan /= num_oscan + + fdata1 = bias_oscans[amp] - mean_oscans[amp] + fmean1 = np.mean(fdata1) + fdata1 -= fmean1 + dmat = np.vstack((reduced_mean_oscan.flatten(), fdata1.flatten())) + covmat = scipy.cov(dmat, rowvar=True) + corr_factor = covmat[0, 1]/covmat[0, 0] + fdiff = fdata1 - corr_factor*reduced_mean_oscan + bias_stats[amp] = BiasStats(np.sqrt(covmat[1, 1]), np.std(fdiff), + corr_factor, + np.mean(bias_oscans[amp])) + #fmean1) + + if make_plots: + f1.suptitle(title) + f2.suptitle(title) + ax1[amp].hist2d(reduced_mean_oscan.flatten(), fdata1.flatten(), + bins=(100, 100), range=((-50, 50), (-50, 50))) + label = 'amp %i, cov/var = %.2f' \ + % (amp, bias_stats[amp].corr_factor) + ax1[amp].text(-40, 40, label, fontsize=6, color='w', + fontweight='bold') + if plot_corr: + ax2[amp].hist(fdiff.flatten(), bins=100, range=(-50, 50), + histtype='step') + else: + ax2[amp].hist(fdata1.flatten(), bins=100, range=(-50, 50), + histtype='step') + + return bias_stats, f1, f2 + + +if __name__ == '__main__': + plt.ion() + bias_files = sorted(glob.glob('/gpfs/slac/lsst/fs1/g/data/jobHarness/jh_archive-test/LCA-11021_RTM/LCA-11021_RTM-002_ETU1/5808D/sflat_raft_acq/v0/38318/S00/ITL-3800C-023-Dev_sflat_bias_*.fits')) + etu1_stats, f1, f2 \ + = correlated_noise(bias_files, make_plots=True, plot_corr=False, + title='Run 5808D, ETU1, S00') + for amp, stats in etu1_stats.items(): + print amp, stats + plt.figure(f1.number) + plt.savefig('ETU1_S00_noise_corr.png') + + bias_files = sorted(glob.glob('/gpfs/slac/lsst/fs1/g/data/jobHarness/jh_archive/LCA-11021_RTM/LCA-11021_RTM-005/6288/sflat_raft_acq/v0/37108/S00/*sflat_bias*.fits')) + rtm005_stats, f1, f2 \ + = correlated_noise(bias_files, make_plots=True, plot_corr=False, + title='Run 6288, RTM-005, S00') + for amp, stats in rtm005_stats.items(): + print amp, stats + plt.figure(f1.number) + plt.savefig('RTM-005_S00_noise_corr.png')