Merge pull request #103 from lsst-camera-dh/LSSTTD-1556_Fe55_gain_sta…

…bility_plots

Lssttd 1556 fe55 gain stability plots

harnessed_jobs/fe55_analysis_BOT/v0/producer_fe55_analysis_BOT.py

@@ -3,9 +3,13 @@
 Producer script for BOT Fe55 analysis.
 """
 import os
+import glob
+from collections import defaultdict
+from bot_eo_analyses import get_analysis_types, run_python_task_or_cl_script
+from fe55_gain_stability import plot_raft_fe55_gains_by_amp, \
+    plot_all_raft_fe55_gains
 from fe55_jh_task import fe55_jh_task
 from gain_stability_jh_task import gain_stability_jh_task
-from bot_eo_analyses import get_analysis_types, run_python_task_or_cl_script
 
 job_dir = os.path.join(os.environ['EOANALYSISJOBSDIR'],
                        'harnessed_jobs', 'fe55_analysis_BOT', 'v0')
@@ -19,3 +23,16 @@
 if 'gainstability' in analysis_types:
     gain_stability_script = os.path.join(job_dir, 'gain_stability_jh_task.py')
     run_python_task_or_cl_script(gain_stability_jh_task, gain_stability_script)
+
+    # Make png files of plots of gains vs MJD for each raft.
+    pickle_files = sorted(glob.glob('*_gain_sequence.pickle'))
+    raft_files = defaultdict(list)
+    for item in pickle_files:
+        raft = os.path.basename(item)[:len('R22')]
+        raft_files[raft].append(item)
+
+    for raft in raft_files:
+        plot_raft_fe55_gains_by_amp(raft_files[raft])
+
+    # Make focal plane summary plot of gain stability, aggregating by CCD.
+    plot_all_raft_fe55_gains(raft_files)

harnessed_jobs/flat_gain_stability_BOT/v0/producer_flat_gain_stability_BOT.py

@@ -3,9 +3,11 @@
 Producer script for BOT flat gain stability analysis.
 """
 import os
+import glob
+from collections import defaultdict
 import siteUtils
 from bot_eo_analyses import get_analysis_types, run_python_task_or_cl_script
-from flat_gain_stability import plot_all_rafts
+from flat_gain_stability import plot_all_rafts, plot_raft_by_amp
 from flat_gain_stability_jh_task import flat_gain_stability_jh_task
 
 if 'gainstability' in get_analysis_types():
@@ -16,4 +18,15 @@
     run_python_task_or_cl_script(flat_gain_stability_jh_task,
                                  flat_gain_stability_script)
 
-    plot_all_rafts(siteUtils.getRunNumber(), y_range=None)
+    # Make png files of plots of gains vs MJD for each raft.
+    pickle_files = sorted(glob.glob('*flat_signal_sequence.pickle'))
+    raft_files = defaultdict(list)
+    for item in pickle_files:
+        raft = os.path.basename(item)[:len('R22')]
+        raft_files[raft].append(item)
+
+    for raft in raft_files:
+        plot_raft_by_amp(raft_files[raft])
+
+    # Make focal plane summary plot of gain stability, aggregating by CCD
+    plot_all_rafts(siteUtils.getRunNumber())

python/bot_eo_analyses.py

@@ -37,6 +37,7 @@
            'append_acq_run',
            'make_title',
            'glob_pattern',
+           'raft_ccd_order',
            'get_mask_files',
            'get_amplifier_gains',
            'medianed_dark_frame',
@@ -182,6 +183,34 @@ def get_mask_files(det_name):
     return mask_files
 
 
+class RaftCcdOrder:
+    """
+    Class to sort eotest output filenames by order of appearance of CCDs
+    in a raft for the standard DM-based focal plane orientation, for
+    science rafts and corner rafts (specifically R44), respectively:
+    S20  S21  S22       SG0
+    S10  S11  S12       SW1
+    S00  S01  S02       SW0  SG1
+    """
+    ccd_order = 'S20 S21 S22 S10 S11 S12 S00 S01 S02 SG0 SW1 SW0 SG1'.split()
+    def _ccd_key(self, filename):
+        """
+        Key function to use in sorted(...).  Filename should be of the
+        form '{raft}_{sensor}_{run}_...'.
+        """
+        ccd_slot = os.path.basename(filename).split('_')[1]
+        return self.ccd_order.index(ccd_slot)
+
+    def sorted(self, file_list):
+        """
+        Return a sorted file list based on CCD slot name.
+        """
+        return sorted(file_list, key=self._ccd_key)
+
+
+raft_ccd_order = RaftCcdOrder()
+
+
 def make_file_prefix(run, component_name):
     """
     Compose the run number and component name into string prefix

python/bot_eo_validators.py

@@ -207,6 +207,11 @@ def validate_fe55(results, det_names):
                 md = dict(DATA_PRODUCT='gain_stability_results')
                 results.append(siteUtils.make_fileref(gain_stability_file,
                                                       metadata=md))
+    # Persist raft-level gain stability plots
+    png_files = glob.glob('*fe55_gain_stability.png')
+    for png_file in png_files:
+        md = dict(DATA_PRODUCT='gain_stability_results')
+        results.append(siteUtils.make_fileref(png_file, metadata=md))
 
     report_missing_data('validate_fe55', missing_det_names)
     if 'gain_stability' in analysis_types:
@@ -713,10 +718,12 @@ def validate_flat_gain_stability(results, det_names):
     report_missing_data('validate_flat_gain_stability', missing_det_names)
 
     unit_id = siteUtils.getUnitId()
-    gain_stability_plot = f'{unit_id}_{run}_flat_gain_stability.png'
-    md = dict(DATA_PRODUCT='flat_gain_stability_plot', LsstId=unit_id)
-    results.append(siteUtils.make_fileref(gain_stability_plot,
-                                          metadata=md))
+    png_files = glob.glob('*flat_gain_stability.png')
+    for png_file in png_files:
+        md = dict(DATA_PRODUCT='flat_gain_stability_plot')
+        if unit_id in png_file:
+            md['LsstId'] = unit_id
+        results.append(siteUtils.make_fileref(png_file, metadata=md))
 
     return results
 

python/fe55_gain_stability.py

@@ -0,0 +1,117 @@
+"""
+Module to make Fe55 gain stability plots.
+"""
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+import siteUtils
+import lsst.eotest.image_utils as imutils
+from bot_eo_analyses import raft_ccd_order, append_acq_run
+
+__all__ = ['plot_raft_fe55_gains_by_amp', 'plot_raft_fe55_gains_by_ccd',
+           'plot_all_raft_fe55_gains']
+
+
+channel = {i: f'C{_}' for i, _ in imutils.channelIds.items()}
+
+
+def plot_raft_fe55_gains_by_amp(raft_files, figsize=(12, 12), y_range=None,
+                                outfile=None):
+    """
+    Plot Fe55 gains for all amps for each CCD in a raft.
+
+    Parameters
+    ----------
+    raft_files: list
+        A list of pickle files containing pandas data frames with
+        the gains for each amp in a CCD as a function of MJD.
+        The filenames are assumed to be of the form
+        '<raft>_<sensor>_<run>_gain_sequence.pickle`, so that the
+        detector name (e.g., 'R22_S11'), raft, and run number can
+        be extracted.
+    figsize: (float, float) [(12, 12)]
+        Size of the figure for each raft.
+    y_range: (float, float) [None]
+        Plotting limts of the y-axis for each plot.  If None, then the
+        y-limits will be scaled by matplotlib to the data.
+    outfile: str [None]
+        Filename of output png file. If None, then a default name of
+        f'{raft}_{run}_fe55_gain_stability.png' will be used.
+    """
+    sorted_raft_files = raft_ccd_order.sorted(raft_files)
+    plt.figure(figsize=figsize)
+    for i, item in enumerate(sorted_raft_files, 1):
+        plt.subplot(3, 3, i)
+        df = pd.read_pickle(item)
+        mjd0 = int(min(df['mjd']))
+        det_name = os.path.basename(item)[:len('R22_S11')]
+        amps = sorted(list(set(df['amp'])))
+        for amp in amps:
+            my_df = df.query(f'amp == {amp}')
+            gains = my_df['gain'].to_numpy()
+            frac = gains/np.mean(gains)
+            plt.scatter(24*(my_df['mjd'] - mjd0), frac, s=2,
+                        label=f'{channel[amp]}')
+        plt.legend(fontsize='x-small', ncol=2)
+        if y_range is not None:
+            plt.ylim(*y_range)
+        plt.title(det_name, fontsize='x-small')
+        plt.xlabel(f'(MJD - {mjd0:d})*24 (hours)')
+        plt.ylabel('gain/mean(gain)')
+    plt.tight_layout(rect=(0, 0, 1, 0.95))
+    tokens = os.path.basename(raft_files[0]).split('_')
+    run = tokens[2]
+    raft = tokens[0]
+    suptitle = append_acq_run(f'Fe55 gain stability, {raft}, Run {run}')
+    plt.suptitle(suptitle)
+    if outfile is None:
+        outfile = f'{raft}_{run}_fe55_gain_stability.png'
+    plt.savefig(outfile)
+
+
+def plot_raft_fe55_gains_by_ccd(raft_files, colors=None, y_range=None):
+    """
+    Plot Fe55 gains for each ccd in a raft, aggregating over amps.
+    """
+    if colors is None:
+        colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+    raft = os.path.basename(raft_files[0])[:len('R22')]
+    mjd0 = None
+    sorted_raft_files = raft_ccd_order.sorted(raft_files)
+    for item, color in zip(sorted_raft_files, colors):
+        sensor = os.path.basename(item)[len('R22_'):len('R22_S11')]
+        df = pd.read_pickle(item)
+        if mjd0 is None:
+            mjd0 = int(min(df['mjd']))
+        x, y = [], []
+        for amp in range(1, 17):
+            my_df = df.query(f'amp == {amp}')
+            time = 24*(my_df['mjd'] - mjd0)
+            gain = my_df['gain']/np.mean(my_df['gain'])
+            x.extend(time)
+            y.extend(gain)
+        plt.scatter(x, y, s=1, color=color, label=sensor)
+    plt.xlabel(f'(MJD - {mjd0})*24 (hours)')
+    plt.ylabel('gain/mean(gain)')
+    plt.legend(fontsize='x-small', ncol=2)
+    plt.title(raft, fontsize='small')
+    if y_range is not None:
+        plt.ylim(*y_range)
+
+
+def plot_all_raft_fe55_gains(raft_files, figsize=(18, 18), y_range=None):
+    """
+    Plot the flat gain stability curve for all 25 rafts in a 5x5 grid.
+    """
+    figure = plt.figure(figsize=figsize)
+    rafts = sorted(list(raft_files.keys()))
+    for i, raft in enumerate(rafts, 1):
+        figure.add_subplot(5, 5, i)
+        plot_raft_fe55_gains_by_ccd(raft_files[raft], y_range=y_range)
+    plt.tight_layout(rect=(0, 0, 1, 0.95))
+    run = os.path.basename(raft_files[rafts[0]][0]).split('_')[2]
+    suptitle = append_acq_run(f'Fe55 gain stability, Run {run}')
+    plt.suptitle(suptitle)
+    unit_id = siteUtils.getUnitId()
+    plt.savefig(f'{unit_id}_{run}_fe55_gain_stability.png')

python/flat_gain_stability.py

@@ -7,19 +7,68 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import pandas as pd
+import lsst.eotest.image_utils as imutils
 import siteUtils
+from bot_eo_analyses import raft_ccd_order, append_acq_run
 
+__all__ = ['plot_raft_by_amp', 'plot_raft', 'plot_all_rafts']
 
-__all__ = ['plot_raft', 'plot_all_rafts']
+
+channel = {i: f'C{_}' for i, _ in imutils.channelIds.items()}
+
+
+def plot_raft_by_amp(raft_files, cut=None, divide_by_flux=True,
+                     figsize=(12, 12), y_range=None, outfile=None):
+    """
+    Plot flat gain stability for each amp individually.
+    """
+    sorted_raft_files = raft_ccd_order.sorted(raft_files)
+    plt.figure(figsize=figsize)
+    for i, item in enumerate(sorted_raft_files, 1):
+        plt.subplot(3, 3, i)
+        df = pd.read_pickle(item)
+        if cut:
+            df = df.query(cut)
+        flux = df['flux']/np.mean(df['flux'])
+        mjd0 = int(min(df['mjd']))
+        det_name = os.path.basename(item)[:len('R22_S11')]
+        for amp in range(1, 17):
+            try:
+                signal = df[f'amp{amp:02d}']
+            except KeyError:
+                continue
+            if divide_by_flux:
+                signal /= flux
+            signal /= np.mean(signal)
+            plt.scatter(24*(df['mjd'] - mjd0), signal, s=2, label=channel[amp])
+        plt.legend(fontsize='x-small', ncol=2)
+        if y_range is not None:
+            plt.ylim(*y_range)
+        plt.title(det_name, fontsize='x-small')
+        plt.xlabel(f'(MJD - {mjd0})*24 (hours)')
+        plt.ylabel('counts/mean(counts)')
+    plt.tight_layout(rect=(0, 0, 1, 0.95))
+    tokens = os.path.basename(raft_files[0]).split('_')
+    run = tokens[2]
+    raft = tokens[0]
+    suptitle = append_acq_run(f'Flat gain stability, {raft}, Run {run}')
+    plt.suptitle(suptitle)
+    if outfile is None:
+        outfile = f'{raft}_{run}_flat_gain_stability.png'
+    plt.savefig(outfile)
 
 
 def plot_raft(raft_files, cut=None, divide_by_flux=True, y_range=None,
               colors=None):
+    """
+    Plot flat gain stability for all ccds in a raft, aggregating over amps.
+    """
     if colors is None:
         colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
     raft = os.path.basename(raft_files[0])[:len('R22')]
     mjd_min = None
-    for item, color in zip(raft_files, colors):
+    sorted_raft_files = raft_ccd_order.sorted(raft_files)
+    for item, color in zip(sorted_raft_files, colors):
         sensor = os.path.basename(item)[len('R22_'):len('R22_S11')]
         raw_df = pd.read_pickle(item)
         df = raw_df.query(cut) if cut is not None else raw_df
@@ -41,7 +90,7 @@ def plot_raft(raft_files, cut=None, divide_by_flux=True, y_range=None,
         plt.scatter(x, y, s=1, color=color, label=sensor)
     plt.xlabel(f'(MJD - {mjd_min})*24 (hours)')
     plt.ylabel('counts/mean(counts)')
-    plt.legend(fontsize='x-small')
+    plt.legend(fontsize='x-small', ncol=2)
     plt.title(raft, fontsize='small')
     if y_range is not None:
         plt.ylim(*y_range)
@@ -55,17 +104,20 @@ def plot_all_rafts(run, results_dir='.', cut=None, y_range=(0.998, 1.002),
     raft_files = defaultdict(list)
     files = sorted(glob.glob(os.path.join(results_dir,
                                           '*flat_signal_sequence.pickle')))
+    if not files:
+        return
     for item in files:
         raft_name = os.path.basename(item)[:len('R22')]
         raft_files[raft_name].append(item)
 
     figure = plt.figure(figsize=figsize)
     rafts = sorted(list(raft_files.keys()))
     for i, raft in enumerate(rafts):
-        ax = figure.add_subplot(5, 5, i+1)
+        figure.add_subplot(5, 5, i+1)
         plot_raft(raft_files[raft], cut=cut, divide_by_flux=divide_by_flux,
                   y_range=y_range)
-    plt.tight_layout()
-    plt.suptitle(f'Run {run}')
+    plt.tight_layout(rect=(0, 0, 1, 0.95))
+    suptitle = append_acq_run(f'Flat gain stability, Run {run}')
+    plt.suptitle(suptitle)
     unit_id = siteUtils.getUnitId()
     plt.savefig(f'{unit_id}_{run}_flat_gain_stability.png')