Fix plotting if cellID data is the same for all trains

1fb92c2a · Karim Ahmed · 7ff7dbae · 1fb92c2a
Commit 1fb92c2a authored 2 years ago by Karim Ahmed
--- a/notebooks/LPD/LPDChar_Darks_NBC.ipynb
+++ b/notebooks/LPD/LPDChar_Darks_NBC.ipynb
@@ -552,7 +552,11 @@
    "            offset = np.nanmedian(data)\n",
    "            noise = np.nanstd(data)\n",
    "            xrange = [np.nanmin(data_g[cap][qm]), np.nanmax(data_g[cap][qm])]\n",
-    "            nbins = int(xrange[1] - xrange[0])\n",
+    "            if len(set(xrange)) == 1:\n",
+    "                xrange = [0, xrange[0]+xrange[0]//2]\n",
+    "                nbins = data_g[cap][qm].shape[0]\n",
+    "            else:\n",
+    "                nbins = int(xrange[1] - xrange[0])\n",
    "\n",
    "            hn, cn = np.histogram(data, bins=nbins, range=xrange)\n",
    "\n",

 %% Cell type:markdown id: tags:
 # LPD Offset, Noise and Dead Pixels Characterization #
 Author: M. Karnevskiy, S. Hauf
 This notebook performs re-characterize of dark images to derive offset, noise and bad-pixel maps. All three types of constants are evaluated per-pixel and per-memory cell.
 The notebook will correctly handle veto settings, but note that if you veto cells you will not be able to use these offsets for runs with different veto settings - vetoed cells will have zero offset.
 The evaluated calibration constants are stored locally and injected in the calibration data base.
 %% Cell type:code id: tags:
 ``` python
 cluster_profile = "noDB" # The ipcluster profile to use
 in_folder = "/gpfs/exfel/exp/FXE/202030/p900121/raw" # path to input data, required
 out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/LPD/" # path to output to, required
 metadata_folder = ""  # Directory containing calibration_metadata.yml when run by xfel-calibrate
 sequence = 0 # sequence files to evaluate
 modules = [-1] # list of modules to evaluate, RANGE ALLOWED
 run_high = 120 # run number in which high gain data was recorded, required
 run_med = 121 # run number in which medium gain data was recorded, required
 run_low = 122 # run number in which low gain data was recorded, required
 karabo_id = "FXE_DET_LPD1M-1" # karabo karabo_id
 karabo_da = ['-1']  # a list of data aggregators names, Default [-1] for selecting all data aggregators
 receiver_id = "{}CH0" # inset for receiver devices
 path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data
 h5path = '/INSTRUMENT/{}/DET/{}:xtdf/image' # path in the HDF5 file to images
 h5path_idx = '/INDEX/{}/DET/{}:xtdf/image' # path in the HDF5 file to images
 use_dir_creation_date = True # use the creation date of the directory for database time derivation
 cal_db_interface = "tcp://max-exfl016:8015#8025" # the database interface to use
 cal_db_timeout = 300000 # timeout on caldb requests"
 local_output = True # output constants locally
 db_output = False # output constants to database
 capacitor_setting = 5 # capacitor_setting for which data was taken
 mem_cells = 512 # number of memory cells used
 bias_voltage = 250 # detector bias voltage
 thresholds_offset_sigma = 3. # bad pixel relative threshold in terms of n sigma offset
 thresholds_offset_hard = [400, 1500] # bad pixel hard threshold
 thresholds_noise_sigma = 7. # bad pixel relative threshold in terms of n sigma noise
 thresholds_noise_hard = [1, 35] # bad pixel hard threshold
 skip_first_ntrains = 10 # Number of first trains to skip
 instrument = "FXE" # instrument name
 ntrains = 100 # number of trains to use
 high_res_badpix_3d = False # plot bad-pixel summary in high resolution
 test_for_normality = False # permorm normality test
 operation_mode = ''  # Detector operation mode, optional
 ```
 %% Cell type:code id: tags:
 ``` python
 import copy
 import os
 import warnings
 from collections import OrderedDict
 from datetime import datetime
 from functools import partial
 warnings.filterwarnings('ignore')
 import dateutil.parser
 import h5py
 import matplotlib
 from ipyparallel import Client
 from IPython.display import Latex, Markdown, display
 matplotlib.use("agg")
 import matplotlib.patches as patches
 import matplotlib.pyplot as plt
 %matplotlib inline
 import numpy as np
 import tabulate
 import yaml
 from iCalibrationDB import Conditions, Constants, Detectors, Versions
 from XFELDetAna.plotting.heatmap import heatmapPlot
 from XFELDetAna.plotting.simpleplot import simplePlot
 from cal_tools.enums import BadPixels
 from cal_tools.plotting import (
    create_constant_overview,
    plot_badpix_3d,
    show_overview,
    show_processed_modules,
 )
 from cal_tools.tools import (
    get_dir_creation_date,
    get_from_db,
    get_notebook_name,
    get_pdu_from_db,
    get_random_db_interface,
    get_report,
    map_gain_stages,
    module_index_to_qm,
    parse_runs,
    run_prop_seq_from_path,
    save_const_to_h5,
    send_to_db,
 )
 ```
 %% Cell type:code id: tags:
 ``` python
 client = Client(profile=cluster_profile)
 view = client[:]
 view.use_dill()
 gains = np.arange(3)
 max_cells = mem_cells
 cells = np.arange(max_cells)
 gain_names = ['High', 'Medium', 'Low']
 if karabo_da[0] == '-1':
    if modules[0] == -1:
        modules = list(range(16))
    karabo_da = ['LPD{:02d}'.format(i) for i in modules]
 else:
    modules = [int(x[-2:]) for x in karabo_da]
 gain_runs = OrderedDict()
 if capacitor_setting == 5:
    gain_runs["high_5pf"] = run_high
    gain_runs["med_5pf"] =  run_med
    gain_runs["low_5pf"] =  run_low
 elif capacitor_setting == 50:
    gain_runs["high_50pf"] = run_high
    gain_runs["med_50pf"] =  run_med
    gain_runs["low_50pf"] =  run_low
 capacitor_settings = [capacitor_setting]
 capacitor_settings = ['{}pf'.format(c) for c in capacitor_settings]
 h5path = h5path.format(karabo_id, receiver_id)
 h5path_idx = h5path_idx.format(karabo_id, receiver_id)
 creation_time = None
 if use_dir_creation_date:
    creation_time = get_dir_creation_date(in_folder, run_high)
    print("Using {} as creation time".format(creation_time))
 run, prop, seq = run_prop_seq_from_path(in_folder)
 cal_db_interface = get_random_db_interface(cal_db_interface)
 display(Markdown('## Evaluated parameters'))
 print('CalDB Interface {}'.format(cal_db_interface))
 print("Proposal: {}".format(prop))
 print("Memory cells: {}/{}".format(mem_cells, max_cells))
 print("Runs: {}, {}, {}".format(run_high, run_med, run_low))
 print("Sequence: {}".format(sequence))
 print("Using DB: {}".format(db_output))
 print("Input: {}".format(in_folder))
 print("Output: {}".format(out_folder))
 print("Bias voltage: {}V".format(bias_voltage))
 ```
 %% Cell type:code id: tags:
 ``` python
 # set everything up filewise
 gmf = map_gain_stages(in_folder, gain_runs, path_template, karabo_da, [sequence])
 gain_mapped_files, total_sequences, total_file_size = gmf
 print(f"Will process a total of {total_sequences} files.")
 ```
 %% Cell type:markdown id: tags:
 ## Data processing
 %% Cell type:code id: tags:
 ``` python
 # the actual characterization
 def characterize_module(cells, bp_thresh, skip_first_ntrains, ntrains, test_for_normality,
                        h5path, h5path_idx, inp):
    import copy
    import h5py
    import numpy as np
    import scipy.stats
    from cal_tools.enums import BadPixels
    def splitOffGainLPD(d):
        msk = np.zeros(d.shape, np.uint16)
        msk[...] = 0b0000111111111111
        data = np.bitwise_and(d, msk)
        msk[...] = 0b0011000000000000
        gain = np.bitwise_and(d, msk)//4096
        gain[gain > 2] = 2
        return data, gain
    filename, channel, gg, cap = inp
    thresholds_offset_hard, thresholds_offset_sigma, thresholds_noise_hard, thresholds_noise_sigma = bp_thresh
    infile = h5py.File(filename, "r")
    h5path = h5path.format(channel)
    h5path_idx = h5path_idx.format(channel)
    count = infile[f"{h5path_idx}/count"][()]
    first = infile[f"{h5path_idx}/first"][()]
    valid = count != 0
    count, first = count[valid], first[valid]
    first_image = int(first[skip_first_ntrains] if first.shape[0] > skip_first_ntrains else 0)
    last_image = int(first_image + np.sum(count[skip_first_ntrains:skip_first_ntrains+ntrains]))
    im = np.array(infile["{}/data".format(h5path, channel)][first_image:last_image, ...])
    cellid = np.squeeze(np.array(infile["{}/cellId".format(h5path, channel)][first_image:last_image, ...]))
    infile.close()
    im, g = splitOffGainLPD(im[:, 0, ...])
    im = im.astype(np.float32)
    im = np.rollaxis(im, 2)
    im = np.rollaxis(im, 2, 1)
    offset = np.zeros((im.shape[0], im.shape[1], cells))
    noise = np.zeros((im.shape[0], im.shape[1], cells))
    normal_test = np.zeros((im.shape[0], im.shape[1], cells))
    for cc in range(cells):
        idx = cellid == cc
        if np.any(idx):
            offset[..., cc] = np.median(im[:, :, idx], axis=2)
            noise[..., cc] = np.std(im[:, :, idx], axis=2)
            if test_for_normality:
                _, normal_test[..., cc] = scipy.stats.normaltest(
                    im[:, :, idx], axis=2)
    # bad pixels
    bp = np.zeros(offset.shape, np.uint32)
    # offset related bad pixels
    offset_mn = np.nanmedian(offset, axis=(0, 1))
    offset_std = np.nanstd(offset, axis=(0, 1))
    bp[(offset < offset_mn-thresholds_offset_sigma*offset_std) |
       (offset > offset_mn+thresholds_offset_sigma*offset_std)] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value
    bp[(offset < thresholds_offset_hard[0]) | (
        offset > thresholds_offset_hard[1])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value
    bp[~np.isfinite(offset)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
    # noise related bad pixels
    noise_mn = np.nanmedian(noise, axis=(0, 1))
    noise_std = np.nanstd(noise, axis=(0, 1))
    bp[(noise < noise_mn-thresholds_noise_sigma*noise_std) |
       (noise > noise_mn+thresholds_noise_sigma*noise_std)] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
    bp[(noise < thresholds_noise_hard[0]) | (
        noise > thresholds_noise_hard[1])] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
    bp[~np.isfinite(noise)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
    idx = cellid == 12
    return offset, noise, channel, gg, cap, bp, im[12, 12, idx], normal_test
 offset_g = OrderedDict()
 noise_g = OrderedDict()
 badpix_g = OrderedDict()
 data_g = OrderedDict()
 ntest_g = OrderedDict()
 gg = 0
 old_cap = None
 start = datetime.now()
 inp = []
 for gain, mapped_files in gain_mapped_files.items():
    cap = gain.split("_")[1]
    if cap != old_cap:
        gg = 0
        old_cap = cap
        offset_g[cap] = OrderedDict()
        noise_g[cap] = OrderedDict()
        badpix_g[cap] = OrderedDict()
        data_g[cap] = OrderedDict()
        ntest_g[cap] = OrderedDict()
    for i in modules:
        qm = module_index_to_qm(i)
        if qm in mapped_files and not mapped_files[qm].empty():
            fname_in = mapped_files[qm].get()
            print("Process file: ", fname_in)
            inp.append((fname_in, i, gg, cap))
    gg+=1
 p = partial(characterize_module, max_cells,
                (thresholds_offset_hard, thresholds_offset_sigma,
                 thresholds_noise_hard, thresholds_noise_sigma),
                skip_first_ntrains, ntrains, test_for_normality,
            h5path, h5path_idx)
 # Don't remove. Used for Debugging.
 #results = list(map(p, inp))
 results = view.map_sync(p, inp)
 for ir, r in enumerate(results):
    offset, noise, i, gg, cap, bp, data, normal = r
    qm = module_index_to_qm(i)
    if qm not in offset_g[cap]:
        offset_g[cap][qm] = np.zeros(
            (offset.shape[0], offset.shape[1], offset.shape[2], 3))
        noise_g[cap][qm] = np.zeros_like(offset_g[cap][qm])
        badpix_g[cap][qm] = np.zeros_like(offset_g[cap][qm], dtype=np.uint32)
        data_g[cap][qm] = np.full((ntrains, 3), np.nan)
        ntest_g[cap][qm] = np.zeros_like(offset_g[cap][qm])
    offset_g[cap][qm][..., gg] = offset
    noise_g[cap][qm][..., gg] = noise
    badpix_g[cap][qm][..., gg] = bp
    data_g[cap][qm][:data.shape[0], gg] = data
    ntest_g[cap][qm][..., gg] = normal
    hn, cn = np.histogram(data, bins=20)
    print(f"{gain_names[gg]} gain, Capacitor {cap}, Module: {qm}. "
          f"Number of processed trains per cell: {data.shape[0]}.")
 ```
 %% Cell type:code id: tags:
 ``` python
 # Read report path and create file location tuple to add with the injection
 proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
 file_loc = 'proposal:{} runs:{} {} {}'.format(proposal, run_low, run_med, run_high)
 report = get_report(metadata_folder)
 ```
 %% Cell type:code id: tags:
 ``` python
 # TODO: add db_module when received from myMDC
 # Create the modules dict of karabo_das and PDUs
 qm_dict = OrderedDict()
 for i, k_da in zip(modules, karabo_da):
    qm = module_index_to_qm(i)
    qm_dict[qm] = {"karabo_da": k_da,
                   "db_module": ""}
 ```
 %% Cell type:code id: tags:
 ``` python
 # Retrieve existing constants for comparison
 clist = ["Offset", "Noise", "BadPixelsDark"]
 old_const = {}
 old_mdata = {}
 dinstance = "LPD1M1"
 detinst = getattr(Detectors, dinstance)
 print('Retrieve pre-existing constants for comparison.')
 for cap in capacitor_settings:
    old_const[cap] = {}
    old_mdata[cap] = {}
    for qm in offset_g[cap].keys():
        old_const[cap][qm] = {}
        old_mdata[cap][qm] = {}
        qm_db = qm_dict[qm]
        karabo_da = qm_db["karabo_da"]
        condition = Conditions.Dark.LPD(memory_cells=max_cells,
                                        bias_voltage=bias_voltage,
                                        capacitor=cap)
        for const in clist:
            constant = getattr(Constants.LPD, const)()
            if not qm_db["db_module"]:
                # This should be used in case of running notebook
                # by a different method other than myMDC which already
                # sends CalCat info.
                qm_db["db_module"] = get_pdu_from_db(karabo_id, [karabo_da], constant,
                                                     condition, cal_db_interface,
                                                     snapshot_at=creation_time)[0]
            data, mdata = get_from_db(karabo_id, karabo_da,
                                      constant,
                                      condition, None,
                                      cal_db_interface,
                                      creation_time=creation_time,
                                      verbosity=2, timeout=cal_db_timeout)
            old_const[cap][qm][const] = data
            if mdata is None or data is None:
                old_mdata[cap][qm][const] = {
                    "timestamp": "Not found",
                    "filepath": None,
                    "h5path": None
                }
            else:
                timestamp = mdata.calibration_constant_version.begin_at.isoformat()
                filepath = os.path.join(
                    mdata.calibration_constant_version.hdf5path,
                    mdata.calibration_constant_version.filename
                )
                h5path = mdata.calibration_constant_version.h5path
                old_mdata[cap][qm][const] = {
                    "timestamp": timestamp,
                    "filepath": filepath,
                    "h5path": h5path
                }
        with open(f"{out_folder}/module_metadata_{qm}.yml","w") as fd:
            yaml.safe_dump(
                {
                    "module": qm,
                    "pdu": qm_db["db_module"],
                    "old-constants": old_mdata[cap][qm]
                }, fd)
 ```
 %% Cell type:code id: tags:
 ``` python
 res = OrderedDict()
 for cap in capacitor_settings:
    res[cap] = OrderedDict()
    for i in modules:
        qm = module_index_to_qm(i)
        res[cap][qm] = {'Offset': offset_g[cap][qm],
                        'Noise': noise_g[cap][qm],
                        'BadPixelsDark': badpix_g[cap][qm]
                        }
 ```
 %% Cell type:code id: tags:
 ``` python
 # Save constants in the calibration DB
 md = None
 for cap in capacitor_settings:
    for qm in res[cap]:
        karabo_da = qm_dict[qm]["karabo_da"]
        db_module = qm_dict[qm]["db_module"]
        # Do not store empty constants
        # In case of 0 trains data_g is initiated with nans and never refilled.
        if np.count_nonzero(~np.isnan(data_g[cap][qm]))==0:
            continue
        for const in res[cap][qm]:
            dconst = getattr(Constants.LPD, const)()
            dconst.data = res[cap][qm][const]
            # set the operating condition
            condition = Conditions.Dark.LPD(memory_cells=max_cells,
                                            bias_voltage=bias_voltage,
                                            capacitor=cap)
            if db_output:
                md = send_to_db(db_module, karabo_id, dconst, condition,
                                file_loc, report_path=report,
                                cal_db_interface=cal_db_interface,
                                creation_time=creation_time,
                                timeout=cal_db_timeout)
            if local_output:
                md = save_const_to_h5(db_module, karabo_id, dconst, condition,
                                      dconst.data, file_loc, report, creation_time, out_folder)
                print(f"Calibration constant {const} is stored locally.\n")
        print("Constants parameter conditions are:\n")
        print(f"• memory_cells: {max_cells}\n• bias_voltage: {bias_voltage}\n"
              f"• capacitor: {cap}\n"
              f"• creation_time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\n")
 ```
 %% Cell type:code id: tags:
 ``` python
 show_processed_modules(
    dinstance=dinstance,
    constants=None,
    mnames=[module_index_to_qm(i) for i in modules],
    mode="position"
 )
 ```
 %% Cell type:markdown id: tags:
 ## Raw pedestal distribution ##
 Distribution of a pedestal (ADUs) over trains for the pixel (12,12), memory cell 12. A median of the distribution is shown in yellow. A standard deviation is shown in red. The green line shows average over all pixels for a given memory cell and gain stage.
 %% Cell type:code id: tags:
 ``` python
 fig, grid = plt.subplots(3, 1, sharex="col", sharey="row", figsize=(10, 7))
 fig.subplots_adjust(wspace=0, hspace=0)
 for cap in capacitor_settings:
    for i in modules:
        qm = module_index_to_qm(i)
        if np.count_nonzero(~np.isnan(data_g[cap][qm])) == 0:
            break
        for gain in range(3):
            data = data_g[cap][qm][:, gain]
            offset = np.nanmedian(data)
            noise = np.nanstd(data)
            xrange = [np.nanmin(data_g[cap][qm]), np.nanmax(data_g[cap][qm])]
-            nbins = int(xrange[1] - xrange[0])
+            if len(set(xrange)) == 1:
+                xrange = [0, xrange[0]+xrange[0]//2]
+                nbins = data_g[cap][qm].shape[0]
+            else:
+                nbins = int(xrange[1] - xrange[0])
            hn, cn = np.histogram(data, bins=nbins, range=xrange)
            grid[gain].hist(data, range=xrange, bins=nbins)
            grid[gain].plot([offset-noise, offset-noise], [0, np.nanmax(hn)],
                            linewidth=1.5, color='red',
                            label='1 $\sigma$ deviation')
            grid[gain].plot([offset+noise, offset+noise],
                            [0, np.nanmax(hn)], linewidth=1.5, color='red')
            grid[gain].plot([offset, offset], [0, 0],
                            linewidth=1.5, color='y', label='median')
            grid[gain].plot([np.nanmedian(offset_g[cap][qm][:, :, 12, gain]),
                             np.nanmedian(offset_g[cap][qm][:, :, 12, gain])],
                            [0, np.nanmax(hn)], linewidth=1.5, color='green',
                            label='average over pixels')
            grid[gain].set_xlim(xrange)
            grid[gain].set_ylim(0, np.nanmax(hn)*1.1)
            grid[gain].set_xlabel("Offset value [ADU]")
            grid[gain].set_ylabel("# of occurance")
            if gain == 0:
                leg = grid[gain].legend(
                    loc='upper center', ncol=3,
                    bbox_to_anchor=(0.1, 0.25, 0.7, 1.0))
            grid[gain].text(820, np.nanmax(hn)*0.4,
                            "{} gain".format(gain_names[gain]), fontsize=20)
            a = plt.axes([.125, .1, 0.775, .8], frame_on=False)
            a.patch.set_alpha(0.05)
            a.set_xlim(xrange)
            plt.plot([offset, offset], [0, 1], linewidth=1.5, color='y')
            plt.xticks([])
            plt.yticks([])
        ypos = 0.9
        x1pos = (np.nanmedian(data_g[cap][qm][:, 0]) +
                 np.nanmedian(data_g[cap][qm][:, 2]))/2.
        x2pos = (np.nanmedian(data_g[cap][qm][:, 2]) +
                 np.nanmedian(data_g[cap][qm][:, 1]))/2.-10
        plt.annotate("", xy=(np.nanmedian(data_g[cap][qm][:, 0]), ypos), xycoords='data',
                     xytext=(np.nanmedian(data_g[cap][qm][:, 2]), ypos), textcoords='data',
                     arrowprops=dict(arrowstyle="<->", connectionstyle="arc3"))
        plt.annotate('{}'.format(np.nanmedian(data_g[cap][qm][:, 0])-np.nanmedian(data_g[cap][qm][:, 2])),
                     xy=(x1pos, ypos), xycoords='data', xytext=(5, 5), textcoords='offset points')
        plt.annotate("", xy=(np.nanmedian(data_g[cap][qm][:, 2]), ypos), xycoords='data',
                     xytext=(np.nanmedian(data_g[cap][qm][:, 1]), ypos), textcoords='data',
                     arrowprops=dict(arrowstyle="<->", connectionstyle="arc3"))
        plt.annotate('{}'.format(np.nanmedian(data_g[cap][qm][:, 2])-np.nanmedian(data_g[cap][qm][:, 1])),
                     xy=(x2pos, ypos), xycoords='data', xytext=(5, 5), textcoords='offset points')
 plt.show()
 ```
 %% Cell type:markdown id: tags:
 ## Normality test ##
 Distributions of raw pedestal values have been tested if they are normally distributed. A normality test have been performed for each pixel and each memory cell. Plots below show histogram of p-Values and a 2D distribution for the  memory cell 12.
 %% Cell type:code id: tags:
 ``` python
 # Loop over capacitor settings, modules, constants
 for cap in capacitor_settings:
    if not test_for_normality:
        print('Normality test was not requested. Flag `test_for_normality` False')
        break
    for i in modules:
        qm = module_index_to_qm(i)
        data = np.copy(ntest_g[cap][qm][:,:,:,:])
        data[badpix_g[cap][qm][:,:,:,:]>0] = 1.01
        hn,cn = np.histogram(data[:,:,:,0], bins=100)
        d = [{'x': np.arange(100)*0.01+0.01,
              'y': np.histogram(data[:,:,:,0], bins=100)[0],
              'drawstyle': 'steps-pre',
              'label' : 'High gain',
              },
             {'x': np.arange(100)*0.01+0.01,
              'y': np.histogram(data[:,:,:,1], bins=100)[0],
              'drawstyle': 'steps-pre',
              'label' : 'Medium gain',
              },
             {'x': np.arange(100)*0.01+0.01,
              'y': np.histogram(data[:,:,:,2], bins=100)[0],
              'drawstyle': 'steps-pre',
              'label' : 'Low gain',
              },
            ]
        fig = plt.figure(figsize=(15,15), tight_layout={'pad': 0.5, 'w_pad': 0.3})
        for gain in range(3):
            ax = fig.add_subplot(221+gain)
            heatmapPlot(data[:,:,12,gain], add_panels=False, cmap='viridis', figsize=(10,10),
                y_label='Rows', x_label='Columns',
                lut_label='p-Value',
                use_axis=ax,
                title='p-Value for cell 12, {} gain'.format(gain_names[gain]) )
        ax = fig.add_subplot(224)
        _ = simplePlot(d, #aspect=1.6,
                              x_label = "p-Value".format(gain),
                              y_label="# of occurance",
                              use_axis=ax,
                               y_log=False, legend='outside-top-ncol3-frame', legend_pad=0.05, legend_size='5%')
        ax.ticklabel_format(style='sci', axis='y', scilimits=(4,6))
 ```
 %% Cell type:raw id: tags:
 .. raw:: latex
    \newpage
 %% Cell type:markdown id: tags:
 ## Single-Cell Overviews ##
 Single cell overviews allow to identify potential effects on all memory cells, e.g. on a sensor level. Additionally, they should serve as a first sanity check on expected behaviour, e.g. if structuring on the ASIC level is visible in the offsets, but otherwise no immediate artifacts are visible.
 %% Cell type:code id: tags:
 ``` python
 cell = 12
 for cap in capacitor_settings:
    for gain in range(3):
        display(
            Markdown('### Cell-12 overview - {} gain'.format(gain_names[gain])))
        fig = plt.figure(figsize=(18, 22) , tight_layout={'pad': 0.1, 'w_pad': 0.1})
        for qm in res[cap]:
            for iconst, const in enumerate(['Offset', 'Noise', 'BadPixelsDark']):
                ax = fig.add_subplot(321+iconst)
                data = res[cap][qm][const][:, :, 12, gain]
                vmax = 1.5 * np.nanmedian(res[cap][qm][const][:, :, 12, gain])
                title = const
                label = '{} value [ADU]'.format(const)
                title = '{} value'.format(const)
                if const == 'BadPixelsDark':
                    vmax = 4
                    bpix_code = data.astype(np.float32)
                    bpix_code[bpix_code == 0] = np.nan
                    title = 'Bad pixel code'
                    label = title
                    cb_labels = ['1 {}'.format(BadPixels.NOISE_OUT_OF_THRESHOLD.name),
                                 '2 {}'.format(BadPixels.OFFSET_NOISE_EVAL_ERROR.name),
                                 '3 {}'.format(BadPixels.OFFSET_OUT_OF_THRESHOLD.name),
                                 '4 {}'.format('MIXED')]
                    heatmapPlot(bpix_code, add_panels=False, cmap='viridis',
                                y_label='Rows', x_label='Columns',
                                lut_label='', vmax=vmax,
                                use_axis=ax, cb_ticklabels=cb_labels, cb_ticks = np.arange(4)+1,
                                title='{}'.format(title))
                    del bpix_code
                else:
                    heatmapPlot(data, add_panels=False, cmap='viridis',
                                y_label='Rows', x_label='Columns',
                                lut_label=label, vmax=vmax,
                                use_axis=ax,
                                title='{}'.format(title))
        for qm in res[cap]:
            for iconst, const in enumerate(['Offset', 'Noise']):
                data = res[cap][qm][const]
                dataBP = np.copy(data)
                dataBP[res[cap][qm]['BadPixelsDark'] > 0] = -1
                x_ranges = [[0, 1500], [0, 40]]
                hn, cn = np.histogram(
                    data[:, :, :, gain], bins=100, range=x_ranges[iconst])
                hnBP, cnBP = np.histogram(dataBP[:, :, :, gain], bins=cn)
                d = [{'x': cn[:-1],
                      'y': hn,
                      'drawstyle': 'steps-pre',
                      'label': 'All data',
                      },
                     {'x': cnBP[:-1],
                      'y': hnBP,
                      'drawstyle': 'steps-pre',
                      'label': 'Bad pixels masked',
                      },
                     ]
                ax = fig.add_subplot(325+iconst)
                _ = simplePlot(d, figsize=(5, 7), aspect=1,
                                    x_label="{} value [ADU]".format(const),
                                    y_label="# of occurance",
                                    title='', legend_pad=0.1, legend_size='10%',
                                    use_axis=ax,
                                    y_log=True, legend='outside-top-2col-frame')
        plt.show()
 ```
 %% Cell type:raw id: tags:
 .. raw:: latex
    \newpage
 %% Cell type:markdown id: tags:
 %% Cell type:code id: tags:
 ``` python
 cols = {BadPixels.NOISE_OUT_OF_THRESHOLD.value: (BadPixels.NOISE_OUT_OF_THRESHOLD.name, '#FF000080'),
        BadPixels.OFFSET_NOISE_EVAL_ERROR.value: (BadPixels.OFFSET_NOISE_EVAL_ERROR.name, '#0000FF80'),
        BadPixels.OFFSET_OUT_OF_THRESHOLD.value: (BadPixels.OFFSET_OUT_OF_THRESHOLD.name, '#00FF0080'),
        BadPixels.OFFSET_OUT_OF_THRESHOLD.value | BadPixels.NOISE_OUT_OF_THRESHOLD.value: ('MIXED', '#DD00DD80')}
 if high_res_badpix_3d:
    display(Markdown("""
    ## Global Bad Pixel Behaviour ##
    The following plots shows the results of a bad pixel evaluation for all evaluated memory cells.
    Cells are stacked in the Z-dimension, while pixels values in x/y are re-binned with a factor of 2.
    This excludes single bad pixels present only in disconnected pixels.
    Hence, any bad pixels spanning at least 4 pixels in the x/y-plane, or across at least two memory cells are indicated.
    Colors encode the bad pixel type, or mixed type.
        """))
    # Switch rebin to 1 for full resolution and
    # no interpolation for badpixel values.
    rebin = 2
    for gain in range(3):
        display(Markdown('### Bad pixel behaviour - {} gain ###'.format(gain_names[gain])))
        for cap in capacitor_settings:
            for mod, data in badpix_g[cap].items():
                plot_badpix_3d(data[...,gain], cols, title='', rebin_fac=rebin)
                ax = plt.gca()
                leg = ax.get_legend()
                leg.set(alpha=0.5)
        plt.show()
 ```
 %% Cell type:raw id: tags:
 .. raw:: latex
    \newpage
 %% Cell type:markdown id: tags:
 ## Summary across tiles ##
 Plots give an overview of calibration constants averaged across tiles. A bad pixel mask is applied. Constants are compared with pre-existing constants retrieved from the calibration database. Differences $\Delta$ between the old and new constants is shown.
 %% Cell type:code id: tags:
 ``` python
 time_summary = []
 for cap, cap_data in old_mdata.items():
    time_summary.append(f"The following pre-existing constants are used for comparison for capacitor setting **{cap}**:")
    for qm, qm_data in cap_data.items():
        time_summary.append(f"- Module {qm}")
        for const, const_data in qm_data.items():
            time_summary.append(f"    - {const} created at {const_data['timestamp']}")
 display(Markdown("\n".join(time_summary)))
 ```
 %% Cell type:code id: tags:
 ``` python
 # Loop over capacitor settings, modules, constants
 for cap in res:
    for qm in res[cap]:
        for gain in range(3):
            display(Markdown('### Summary across tiles - {} gain'.format(gain_names[gain])))
            for const in res[cap][qm]:
                data = np.copy(res[cap][qm][const][:, :, :, gain])
                label = 'Fraction of bad pixels'
                if const != 'BadPixelsDark':
                    data[badpix_g[cap][qm][:, :, :, gain] > 0] = np.nan
                    label = '{} value [ADU]'.format(const)
                else:
                    data[data>0] = 1.0
                data = data.reshape(
                    int(data.shape[0] / 32),
                    32,
                    int(data.shape[1] / 128),
                    128,
                    data.shape[2])
                data = np.nanmean(data, axis=(1, 3)).swapaxes(
                    0, 2).reshape(512, 16)
                fig = plt.figure(figsize=(15, 6))
                ax = fig.add_subplot(121)
                _ = heatmapPlot(data[:510, :], add_panels=True,
                                y_label='Momery Cell ID', x_label='Tile ID',
                                lut_label=label, use_axis=ax,
                                panel_y_label=label, panel_x_label=label,
                                cmap='viridis',  # cb_loc='right',cb_aspect=15,
                                x_ticklabels=np.arange(16)+1,
                                x_ticks=np.arange(16)+0.5)
                if old_const[cap][qm][const] is not None:
                    ax = fig.add_subplot(122)
                    dataold = np.copy(old_const[cap][qm][const][:, :, :, gain])
                    label = '$\Delta$ {}'.format(label)
                    if const != 'BadPixelsDark':
                        if old_const[cap][qm]['BadPixelsDark'] is not None:
                            dataold[old_const[cap][qm]['BadPixelsDark'][:, :, :, gain] > 0] = np.nan
                        else:
                            dataold[:] = np.nan
                    else:
                        dataold[dataold>0]=1.0
                    dataold = dataold.reshape(
                        int(dataold.shape[0] / 32),
                        32,
                        int(dataold.shape[1] / 128),
                        128,
                        dataold.shape[2])
                    dataold = np.nanmean(dataold, axis=(
                        1, 3)).swapaxes(0, 2).reshape(512, 16)
                    dataold = dataold - data
                    _ = heatmapPlot(dataold[:510, :], add_panels=True,
                                    y_label='Momery Cell ID', x_label='Tile ID',
                                    lut_label=label, use_axis=ax,
                                    panel_y_label=label, panel_x_label=label,
                                    cmap='viridis',  # cb_loc='right',cb_aspect=15,
                                    x_ticklabels=np.arange(16)+1,
                                    x_ticks=np.arange(16)+0.5)
            plt.show()
 ```
 %% Cell type:raw id: tags:
 .. raw:: latex
    \newpage
 %% Cell type:markdown id: tags:
 ## Variation of offset and noise across Tiles and ASICs ##
 The following plots show a standard deviation $\sigma$ of the calibration constant. The plot of standard deviations across tiles show pixels of one tile ($128 \times 32$). Value for each pixel shows a standard deviation across 16 tiles. The standard deviation across ASICs are shown overall tiles. The plot shows pixels of one ASIC ($16 \times 32$), where the value shows a standard deviation across all ACIS of the module.
 %% Cell type:code id: tags:
 ``` python
 # Loop over capacitor settings, modules, constants
 for cap in res:
    for qm in res[cap]:
        for gain in range(3):
            display(Markdown('### Variation of offset and noise across ASICs - {} gain'.format(gain_names[gain])))
            fig = plt.figure(figsize=(15, 6))
            for iconst, const in enumerate(['Offset', 'Noise']):
                data = np.copy(res[cap][qm][const][:, :, :, gain])
                data[badpix_g[cap][qm][:, :, :, gain] > 0] = np.nan
                label = '$\sigma$ {} [ADU]'.format(const)
                dataA = np.nanmean(data, axis=2)  # average over cells
                dataA = dataA.reshape(8, 32, 16, 16)
                dataA = np.nanstd(dataA, axis=(0, 2))  # average across ASICs
                ax = fig.add_subplot(121+iconst)
                _ = heatmapPlot(dataA, add_panels=True,
                                y_label='rows', x_label='columns',
                                lut_label=label, use_axis=ax,
                                panel_y_label=label, panel_x_label=label,
                                cmap='viridis'
                                )
            plt.show()
 ```
 %% Cell type:code id: tags:
 ``` python
 # Loop over capacitor settings, modules, constants
 for cap in res:
    for qm in res[cap]:
        for gain in range(3):
            display(Markdown('### Variation of offset and noise across tiles - {} gain'.format(gain_names[gain])))
            fig = plt.figure(figsize=(15, 6))
            for iconst, const in enumerate(['Offset', 'Noise']):
                data = np.copy(res[cap][qm][const][:, :, :, gain])
                data[badpix_g[cap][qm][:, :, :, gain] > 0] = np.nan
                label = '$\sigma$ {} [ADU]'.format(const)
                dataT = data.reshape(
                    int(data.shape[0] / 32),
                    32,
                    int(data.shape[1] / 128),
                    128,
                    data.shape[2])
                dataT = np.nanstd(dataT, axis=(0, 2))
                dataT = np.nanmean(dataT, axis=2)
                ax = fig.add_subplot(121+iconst)
                _ = heatmapPlot(dataT, add_panels=True,
                                y_label='rows', x_label='columns',
                                lut_label=label, use_axis=ax,
                                panel_y_label=label, panel_x_label=label,
                                cmap='viridis')
            plt.show()
 ```
 %% Cell type:raw id: tags:
 .. raw:: latex
    \newpage
 %% Cell type:markdown id: tags:
 ## Aggregate values and per cell behaviour ##
 The following tables and plots give an overview of statistical aggregates for each constant, as well as per-cell behavior, averaged across pixels.
 %% Cell type:code id: tags:
 ``` python
 # Loop over capacitor settings, modules, constants
 for cap in res:
    for qm in res[cap]:
        for gain in range(3):
            display(Markdown('### Mean over pixels - {} gain'.format(gain_names[gain])))
            fig = plt.figure(figsize=(9,11))
            for iconst, const in enumerate(res[cap][qm]):
                ax = fig.add_subplot(311+iconst)
                data = res[cap][qm][const][:,:,:510,gain]
                if const == 'BadPixelsDark':
                    data[data>0] = 1.0
                dataBP = np.copy(data)
                dataBP[badpix_g[cap][qm][:,:,:510,gain]>0] = -10
                data = np.nanmean(data, axis=(0,1))
                dataBP = np.nanmean(dataBP, axis=(0,1))
                d = [{'y': data,
                      'x': np.arange(data.shape[0]),
                      'drawstyle': 'steps-mid',
                      'label' : 'All data'
                     }
                    ]
                if const != 'BadPixelsDark':
                    d.append({'y': dataBP,
                      'x': np.arange(data.shape[0]),
                      'drawstyle': 'steps-mid',
                      'label' : 'good pixels only'
                     })
                    y_title = "{} value [ADU]".format(const)
                    title = "{} value, {} gain".format(const, gain_names[gain])
                else:
                    y_title = "Fraction of Bad Pixels"
                    title = "Fraction of Bad Pixels, {} gain".format(gain_names[gain])
                data_min = np.min([data, dataBP])if const != 'BadPixelsDark' else np.min([data])
                data_max = np.max([data[20:], dataBP[20:]])
                data_dif = data_max - data_min
                local_max = np.max([data[200:300], dataBP[200:300]])
                frac = 0.35
                new_max = (local_max - data_min*(1-frac))/frac
                new_max = np.max([data_max, new_max])
                _ = simplePlot(d, figsize=(10,10), aspect=2, xrange=(-12, 510),
                                  x_label = 'Memory Cell ID',
                                  y_label=y_title, use_axis=ax,
                                  title=title,
                                  title_position=[0.5, 1.15],
                                  inset='xy-coord-right', inset_x_range=(0,20), inset_indicated=True,
                                  inset_labeled=True, inset_coord=[0.2,0.5,0.6,0.95],
                                    inset_lw = 1.0, y_range = [data_min-data_dif*0.05, new_max+data_dif*0.05],
                                  y_log=False, legend='outside-top-ncol2-frame', legend_size='18%',
                                     legend_pad=0.00)
                plt.tight_layout(pad=1.08, h_pad=0.35)
            plt.show()
 ```
 %% Cell type:raw id: tags:
 .. raw:: latex
    \newpage
 %% Cell type:markdown id: tags:
 ## Summary tables ##
 The following tables show summary information for the evaluated module. Values for currently evaluated constants are compared with values for pre-existing constants retrieved from the calibration database.
 %% Cell type:code id: tags:
 ``` python
 table = []
 bits = [BadPixels.NOISE_OUT_OF_THRESHOLD, BadPixels.OFFSET_OUT_OF_THRESHOLD, BadPixels.OFFSET_NOISE_EVAL_ERROR]
 for cap in res:
    for qm in res[cap]:
        for gain in range(3):
            l_data = []
            l_data_old = []
            data = np.copy(res[cap][qm]['BadPixelsDark'][:,:,:,gain])
            l_data.append(len(data[data>0].flatten()))
            for bit in bits:
                l_data.append(np.count_nonzero(badpix_g[cap][qm][:,:,:,gain] & bit.value))
            if old_const[cap][qm]['BadPixelsDark'] is not None:
                old_const[cap][qm]['BadPixelsDark'] = old_const[cap][qm]['BadPixelsDark'].astype(np.uint32)
                dataold = np.copy(old_const[cap][qm]['BadPixelsDark'][:, :, :, gain])
                l_data_old.append(len(dataold[dataold>0].flatten()))
                for bit in bits:
                    l_data_old.append(np.count_nonzero(old_const[cap][qm]['BadPixelsDark'][:, :, :, gain] & bit.value))
            l_data_name = ['All bad pixels', 'NOISE_OUT_OF_THRESHOLD',
                           'OFFSET_OUT_OF_THRESHOLD', 'OFFSET_NOISE_EVAL_ERROR']
            l_threshold = ['', f'{thresholds_noise_sigma}', f'{thresholds_offset_sigma}',
                           f'{thresholds_offset_hard}/{thresholds_noise_hard}']
            for i in range(len(l_data)):
                line = [f'{l_data_name[i]}, gain {gain_names[gain]}', l_threshold[i], l_data[i]]
                if old_const[cap][qm]['BadPixelsDark'] is not None:
                    line += [l_data_old[i]]
                else:
                    line += ['-']
                table.append(line)
            table.append(['', '', '', ''])
 display(Markdown('''
 ### Number of bad pixels ###
 One pixel can be bad for different reasons, therefore, the sum of all types of bad pixels can be more than the number of all bad pixels.
 '''))
 if len(table)>0:
    md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
                                     headers=["Pixel type", "Threshold",
                                              "New constant", "Old constant"])))
 ```
 %% Cell type:code id: tags:
 ``` python
 header = ['Parameter',
          "New constant", "Old constant ",
          "New constant", "Old constant ",
          "New constant", "Old constant "]
 for const in ['Offset', 'Noise']:
    table = [['','High gain', 'High gain', 'Medium gain', 'Medium gain', 'Low gain', 'Low gain']]
    for cap in res:
        for qm in res[cap]:
            data = np.copy(res[cap][qm][const])
            data[res[cap][qm]['BadPixelsDark']>0] = np.nan
            if old_const[cap][qm][const] is not None and old_const[cap][qm]['BadPixelsDark'] is not None :
                dataold = np.copy(old_const[cap][qm][const])
                dataold[old_const[cap][qm]['BadPixelsDark']>0] = np.nan
            f_list = [np.nanmedian, np.nanmean, np.nanstd, np.nanmin, np.nanmax]
            n_list = ['Median', 'Mean', 'Std', 'Min', 'Max']
            for i, f in enumerate(f_list):
                line = [n_list[i]]
                for gain in range(3):
                    line.append('{:6.1f}'.format(f(data[...,gain])))
                    if old_const[cap][qm][const] is not None and old_const[cap][qm]['BadPixelsDark'] is not None:
                        line.append('{:6.1f}'.format(f(dataold[...,gain])))
                    else:
                        line.append('-')
                table.append(line)
    display(Markdown('### {} [ADU], good pixels only ###'.format(const)))
    md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=header)))
 ```