Fixed BadPixelsCS swaped axes.

e5b7489c · Ivana Klackova · 7d2fb9ce · e5b7489c
Commit e5b7489c authored 1 year ago by Ivana Klackova
--- a/notebooks/AGIPD/CS_Characterization_unequalClockStep_Summary.ipynb
+++ b/notebooks/AGIPD/CS_Characterization_unequalClockStep_Summary.ipynb
@@ -300,7 +300,7 @@
    "            if const == \"SlopesCS\":\n",
    "                dbconst.data = slope_dict_to_arr(sanitised_const[mod])\n",
    "            else:\n",
-    "                dbconst.data = BPmap[mod]\n",
+    "                dbconst.data = BPmap[mod].swapaxes(1,2) \n",
    "\n",
    "            md = send_to_db(pdu, karabo_id, dbconst, condition,\n",
    "                            file_loc, report, cal_db_interface,\n",

 %% Cell type:markdown id: tags:

 # Current Source Characterisation Summary

 This notebook is used as a dependency notebook for current source characterisation to provide summary for all modules of the AGIPD.

 %% Cell type:code id: tags:

 ``` python
 in_folder = "" # in this notebook, in_folder is not used as the data source is in the destination folder
 out_folder = "/gpfs/exfel/exp/SPB/202330/p900340/scratch/CS_Processing/test/230914/"  # the folder to output to, required
 metadata_folder = ""  # Directory containing calibration_metadata.yml when run by xfel-calibrate
 proc_folder = "" # Path to corrected image data used to create histograms and validation plots
 raw_folder = '/gpfs/exfel/exp/SPB/202330/p900340/raw/'  # folder of raw data. This is used to save information of source data of generated constants, required
 dark_run = 5

 karabo_id_control = "SPB_IRU_AGIPD1M1"  # karabo-id for the control device e.g. "MID_EXP_AGIPD1M1", or "SPB_IRU_AGIPD1M1"
 karabo_id = "SPB_DET_AGIPD1M-1"
 ctrl_source_template = '{}/MDL/FPGA_COMP' # path to control information

 creation_time = "" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC e.g. "2022-06-28 13:00:00"
 creation_date_offset = "00:00:00" # add an offset to creation date, e.g. to get different constants
 cal_db_interface = "tcp://max-exfl-cal002:8015#8045"
 db_output = False

 # Detector conditions
 bias_voltage = -1 # detector bias voltage, use -1 to use value stored in slow data.
 mem_cells = -1  # number of memory cells used, use -1 to use value stored in slow data.
 acq_rate = -1. # the detector acquisition rate, use -1. to use value stored in slow data.
 gain_setting = -1 # the gain setting, use -1 to use value stored in slow data.
 ```

 %% Cell type:code id: tags:

 ``` python
 import warnings

 import h5py
 import matplotlib.pyplot as plt
 import matplotlib.gridspec as gridspec
 import matplotlib.cm as cm
 import numpy as np

 import tabulate
 from datetime import timedelta
 from dateutil import parser
 import os

 from cal_tools.agipdlib import AgipdCtrl

 from cal_tools.ana_tools import get_range
 from cal_tools.enums import BadPixels
 from cal_tools.tools import (
    module_index_to_qm,
    calcat_creation_time,
    get_report,
    get_pdu_from_db,
    send_to_db
 )

 from iCalibrationDB import Conditions, Constants

 from extra_data import RunDirectory
 from extra_geom import AGIPD_1MGeometry, AGIPD_500K2GGeometry

 from IPython.display import Latex, display
 from XFELDetAna.plotting.simpleplot import simplePlot

 %matplotlib inline
 warnings.filterwarnings('ignore')
 ```

 %% Cell type:code id: tags:

 ``` python
 # Get operation conditions
 ctrl_source = ctrl_source_template.format(karabo_id_control)
 run_folder = f'{raw_folder}/r{dark_run:04d}/'
 raw_dc = RunDirectory(run_folder)

 # Read operating conditions from AGIPD00 files
 instrument_src_mod = [
    s for s in list(raw_dc.all_sources) if "0CH" in s][0]
 ctrl_src = [
    s for s in list(raw_dc.all_sources) if ctrl_source in s][0]
 ```

 %% Cell type:code id: tags:

 ``` python
 # Evaluate creation time
 creation_time = calcat_creation_time(raw_folder, dark_run, creation_time)
 offset = parser.parse(creation_date_offset)
 delta = timedelta(hours=offset.hour, minutes=offset.minute, seconds=offset.second)
 creation_time += delta
 print(f"Creation time: {creation_time}")

 agipd_cond = AgipdCtrl(
    run_dc=raw_dc,
    image_src=instrument_src_mod,
    ctrl_src=ctrl_src,
    raise_error=False,  # to be able to process very old data without mosetting value
 )
 if mem_cells == -1:
    mem_cells = agipd_cond.get_num_cells()
 if mem_cells is None:
    raise ValueError(f"No raw images found in {run_folder}")
 if acq_rate == -1.:
    acq_rate = agipd_cond.get_acq_rate()
 if gain_setting == -1:
    gain_setting = agipd_cond.get_gain_setting(creation_time)
 if bias_voltage == -1:
    bias_voltage = agipd_cond.get_bias_voltage(karabo_id_control)

 # Evaluate detector instance for mapping
 instrument = karabo_id.split("_")[0]
 nmods = 8 if instrument == 'HED' else 16

 print(f"Using {creation_time} as creation time")
 print(f"Operating conditions are:\n• Bias voltage: {bias_voltage}\n• Memory cells: {mem_cells}\n"
      f"• Acquisition rate: {acq_rate}\n• Gain setting: {gain_setting}"
 )
 ```

 %% Cell type:code id: tags:

 ``` python
 # Sanitised constants keys:
 # mH: high gain slope
 # mM: medium gain slope
 # mL: low gain slope
 #
 # bH: high gain intercept
 # bM: medium gain intercept
 # bL: low gain intercept
 #
 # H-M: ratio of high gain and medium gain slope
 # M-L: ratio of medium gain and low gain slope

 # Consts saved to DB as (8,mem_cells,128,512)
 # order: ['mH', 'mM', 'mL', 'bH', 'bM', 'bL']

 keys_fit = ['m', 'b']
 gs = ['H', 'M', 'L']
 keys_ratio = ['H-M', 'M-L']

 labels = {'m': 'Slope (m)',
          'b': 'Intercept (b)',
          'H-M': "High-to-Medium Ratio",
          'M-L': 'Medium-to-Low Ratio',
          'mask': 'Fraction of bad pixels over cells' }

 fit_data = {}
 ratios = {}
 BPmap = {}
 sanitised_const = {}

 modules = []
 karabo_da = []
 for mod in range(nmods):
    qm = module_index_to_qm(mod)
    ratios[mod] = {}
    fit_data[mod] = {}
    sanitised_const[mod] = {}
    constants_file = f'{out_folder}/CSconst_{karabo_id}_M{mod}.h5'

    if os.path.exists(constants_file):
        print(f'Data available for module {qm}')

        with h5py.File(constants_file, 'r') as hf:
            BPmap[mod] = hf['/BadPixels/data'][()].swapaxes(1,2)

            for key in keys_fit:
                fit_data[mod][key] = {}
                for g in range(0,3):
                    # loop 1st through keys, then through gains
                    # to have right order of constants for DB injection
                    fit_data[mod][key][g] = hf[f'/SanitizedConsts/{g}/{key}/data'][()].swapaxes(1,2)
                    sanitised_const[mod][f'{key}{gs[g]}'] = hf[f'/SanitizedConsts/{g}/{key}/data'][()]

            for key in keys_ratio:
                ratios[mod][key]= hf[f'/SanitizedConsts/Ratios/{key}/data'][()].swapaxes(1,2)
                sanitised_const[mod][key] = hf[f'/SanitizedConsts/Ratios/{key}/data'][()]

        modules.append(mod)
        karabo_da.append(f"AGIPD{mod:02d}")
    else:
        print(f"No fit data available for module {qm}")

 # Change the order of dict keys to maintain compatibility for the rest of code
 fit_data = {mod: {g: {f: fit_data[mod][f][g] for f in keys_fit
                     }
                  for g in range(0,3)
                 }
            for mod in modules
           }
 ```

 %% Cell type:code id: tags:

 ``` python
 def slope_dict_to_arr(d):
    """Convert dictionary to numpy array."""

    arr = np.zeros((8,mem_cells,128,512), np.float32)
    for i, key in enumerate(d):
        arr[i,...] = d[key]

    return arr
 ```

 %% Cell type:code id: tags:

 ``` python
 proposal = list(filter(None, raw_folder.strip('/').split('/')))[-2]
 file_loc = f'Proposal: {proposal}, Run: {dark_run}'

 report = get_report(metadata_folder)
 ```

 %% Cell type:markdown id: tags:

 ## Injection to DB

 %% Cell type:code id: tags:

 ``` python
 if not db_output:
    print('Injection to DB not requested.')
 ```

 %% Cell type:code id: tags:

 ``` python
 md = None

 # set the operating condition
 condition = Conditions.Dark.AGIPD(memory_cells=mem_cells,
                                  bias_voltage=bias_voltage,
                                  acquisition_rate=acq_rate,
                                  gain_setting=gain_setting
                                 )

 db_modules = get_pdu_from_db(karabo_id, karabo_da, Constants.AGIPD.SlopesCS(),
                             condition, cal_db_interface,
                             snapshot_at=creation_time)

 cond_dev = {"Memory cells": [352, 352],
            "Acquisition rate": [4.5, 4.5]
           }

 if db_output:
    for parm in condition.parameters:
        if parm.name in cond_dev:
            parm.lower_deviation = cond_dev[parm.name][0]
            parm.upper_deviation = cond_dev[parm.name][1]
            print(f'Condition for {parm.name} deviation updated to {cond_dev[parm.name]}.')

    for mod, pdu in zip(modules, db_modules):
        for const in ["SlopesCS", "BadPixelsCS"]:
            dbconst = getattr(Constants.AGIPD, const)()

            if const == "SlopesCS":
                dbconst.data = slope_dict_to_arr(sanitised_const[mod])
            else:
-                dbconst.data = BPmap[mod]
+                dbconst.data = BPmap[mod].swapaxes(1,2)

            md = send_to_db(pdu, karabo_id, dbconst, condition,
                            file_loc, report, cal_db_interface,
                            creation_time=creation_time)

    print("Constants parameter conditions are:\n")
    print(f"• memory_cells: {np.diff(cond_dev['Memory cells'])[0]}-{cond_dev['Memory cells'][0]}\n"
          f"• bias_voltage: {bias_voltage}\n"
          f"• acquisition_rate: {np.diff(cond_dev['Acquisition rate'])[0]}-{cond_dev['Acquisition rate'][0]}\n"
          f"• gain_setting: {gain_setting}\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
 #Define AGIPD geometry
 if instrument == "HED":
    geom = AGIPD_500K2GGeometry.from_origin()
 else:
    geom = AGIPD_1MGeometry.from_quad_positions(quad_pos=[
        (-525, 625),
        (-550, -10),
        (520, -160),
        (542.5, 475),
    ])
 ```

 %% Cell type:code id: tags:

 ``` python
 # Create the arrays that will be used for figures.
 # Each array correponds to the data for all processed modules.

 pixel_range = [0,0,512,128]

 # const_data contains ratios of slopes and BP
 const_data = {}
 for key in keys_ratio:
    const_data[key] = np.full((nmods, mem_cells, 512, 128), np.nan)
    for mod in modules:
        if key in ratios[mod]:
            const_data[key][mod,:,pixel_range[0]:pixel_range[2],
                               pixel_range[1]:pixel_range[3]] = ratios[mod][key]

 const_data['mask'] = np.full((nmods, mem_cells, 512, 128), np.nan)
 for mod in modules:
    const_data['mask'][mod,:,pixel_range[0]:pixel_range[2],
                    pixel_range[1]:pixel_range[3]] = BPmap[mod]

 # fit_const_data contains fitted slopes and intercepts for each gain stage
 fit_const_data = {}
 for g in range(0,3):
    fit_const_data[g] = {}
    for key in keys_fit:
        fit_const_data[g][key] = np.full((nmods, mem_cells, 512, 128), np.nan)
        for mod in modules:
            if key in fit_data[mod][g].keys():
                fit_const_data[g][key][mod,:,pixel_range[0]:pixel_range[2],
                               pixel_range[1]:pixel_range[3]] = fit_data[mod][g][key]
 ```

 %% Cell type:markdown id: tags:

 ## Summary across pixels ##

 %% Cell type:code id: tags:

 ``` python
 g_label = ['High gain', 'Medium gain', 'Low gain']
 for g in range(0,3):
    gs = gridspec.GridSpec(1, 2)
    fig = plt.figure(figsize=(15, 7))
    plt.suptitle(f'{g_label[g]}', fontsize=16)

    for i, key in enumerate(fit_const_data[g].keys()):
        data = np.nanmean(fit_const_data[g][key], axis=1)
        vmin, vmax = get_range(data, 5)
        ax = fig.add_subplot(gs[0, i])
        geom.plot_data_fast(data,
                                vmin=vmin, vmax=vmax, ax=ax, cmap="jet", figsize=(13,7),
                                colorbar={'shrink': 1,
                                          'pad': 0.04,
                                          'fraction': 0.1
                                         })
        colorbar = ax.images[0].colorbar
        ax.set_title(labels[key], fontsize=14)
        ax.set_xlabel('Columns', fontsize=12)
        ax.set_ylabel('Rows', fontsize=12)

    plt.show()
 ```

 %% Cell type:code id: tags:

 ``` python
 gs = gridspec.GridSpec(2, 2)
 fig = plt.figure(figsize=(15, 15), tight_layout=True)

 for i, key in enumerate(const_data.keys()):
    if key=='mask':
        ax = fig.add_subplot(gs[0, :])
        data = np.nanmean(const_data[key]>0, axis=1)
        vmin, vmax = (0,1)
    else:
        ax = fig.add_subplot(gs[1, i])
        data = np.nanmean(const_data[key], axis=1)
        vmin, vmax = get_range(data, 5)
    geom.plot_data_fast(data,
                            vmin=vmin, vmax=vmax, ax=ax, cmap="jet", figsize=(13,7),
                            colorbar={'shrink': 1,
                                      'pad': 0.04,
                                      'fraction': 0.1
                                     })
    colorbar = ax.images[0].colorbar
    ax.set_title(labels[key], fontsize=14)
    ax.set_xlabel('Columns', fontsize=12)
    ax.set_ylabel('Rows', fontsize=12)

 plt.show()
 ```

 %% Cell type:markdown id: tags:

 ## Summary across cells ##

 Good pixels only.

 %% Cell type:code id: tags:

 ``` python
 for key in const_data.keys():
    data = np.copy(const_data[key])
    if key=='mask':
        data = data>0
    else:
        data[const_data['mask']>0] = np.nan

    d = []
    for i in range(nmods):
        d.append({'x': np.arange(data[i].shape[0]),
                  'y': np.nanmean(data[i], axis=(1,2)),
                  'drawstyle': 'steps-pre',
                  'label': f'{i}',
                  'linewidth': 2,
                  'linestyle': '--' if i>7 else '-'
                  })

    fig = plt.figure(figsize=(15, 6))
    ax = fig.add_subplot(111)

    _ = simplePlot(d, xrange=(-12, 510),
                        x_label='Memory Cell ID',
                        y_label=labels[key],
                        use_axis=ax,
                        legend='top-left-frame-ncol8',)
    ylim = ax.get_ylim()
    ax.set_ylim(ylim[0], ylim[1] + np.abs(ylim[1]-ylim[0])*0.2)
    ax.grid()
 ```

 %% Cell type:code id: tags:

 ``` python
 table = []
 for mod in modules:

    table.append((mod,
                  f"{np.nanmean(ratios[mod]['H-M']):0.1f} +- {np.nanstd(ratios[mod]['H-M']):0.2f}",
                  f"{np.nanmean(ratios[mod]['M-L']):0.1f} +- {np.nanstd(ratios[mod]['M-L']):0.2f}",
                  f"{np.nanmean(BPmap[mod]>0)*100:0.1f} ({np.nansum(BPmap[mod]>0)})"
                        ))

 all_HM = []
 all_ML = []
 for mod in modules:
    all_HM.extend(ratios[mod]['H-M'])
    all_ML.extend(ratios[mod]['M-L'])

 all_HM = np.array(all_HM)
 all_ML = np.array(all_ML)

 all_MSK = np.array([list(msk) for msk in BPmap.values()])

 table.append(('overall',
              f"{np.nanmean(all_HM):0.1f} +- {np.nanstd(all_HM):0.2f}",
              f"{np.nanmean(all_ML):0.1f} +- {np.nanstd(all_ML):0.2f}",
              f"{np.nanmean(all_MSK>0)*100:0.1f} ({np.nansum(all_MSK>0)})"
                    ))

 md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
                                     headers=["Module", "Ratio High-Medium", "Ratio Medium-Low","Bad pixels [%(Count)]"])))
 ```

 %% Cell type:markdown id: tags:

 ## Summary gain ratios (good pixels only) + histograms

 %% Cell type:code id: tags:

 ``` python
 gs = gridspec.GridSpec(2, 2)
 fig = plt.figure(figsize=(15, 15))

 colors = cm.rainbow(np.linspace(0, 1, nmods))
 for i, key in enumerate(const_data.keys()):
    if key != 'mask':
        const_data[key][const_data['mask']>0] = np.nan
        data = np.nanmean(const_data[key], axis=1)
        vmin, vmax = get_range(data, 5)
        ax = fig.add_subplot(gs[0, i])
        geom.plot_data_fast(data,
                                vmin=vmin, vmax=vmax, ax=ax, cmap="jet", figsize=(13,7),
                                colorbar={'shrink': 1,
                                          'pad': 0.04,
                                          'fraction': 0.1
                                         })
        colorbar = ax.images[0].colorbar
        ax.set_title(labels[key], fontsize=14)
        ax.set_xlabel('Columns', fontsize=12)
        ax.set_ylabel('Rows', fontsize=12)

        ax = fig.add_subplot(gs[1,i])
        for nmod in range(const_data[key].shape[0]):
            h, e = np.histogram(const_data[key][nmod,...].flatten(), bins=100, range=(vmin, vmax))
            ax.plot(e[:-1], h, color=colors[nmod],linewidth=2, label=f'{nmod}', alpha=0.8)
            ax.set_xlabel(f'{labels[key]}')
            ax.set_ylabel('Counts')
        ax.grid()
        ax.legend()



 plt.show()
 ```