{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# LPD Offset, Noise and Dead Pixels Characterization #\n", "\n", "Author: M. Karnevskiy, S. Hauf\n", "\n", "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.\n", "\n", "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.\n", "\n", "The evaluated calibration constants are stored locally and injected in the calibration data base.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "in_folder = \"/gpfs/exfel/exp/FXE/202304/p003338/raw\" # path to input data, required\n", "out_folder = \"/gpfs/exfel/data/scratch/kluyvert/lpd-darks-p3338-r133-134-135/\" # path to output to, required\n", "metadata_folder = \"\" # Directory containing calibration_metadata.yml when run by xfel-calibrate\n", "modules = [-1] # list of modules to evaluate, RANGE ALLOWED\n", "run_high = 133 # run number in which high gain data was recorded, required\n", "run_med = 134 # run number in which medium gain data was recorded, required\n", "run_low = 135 # run number in which low gain data was recorded, required\n", "\n", "karabo_id = \"FXE_DET_LPD1M-1\" # karabo karabo_id\n", "karabo_da = ['-1'] # a list of data aggregators names, Default [-1] for selecting all data aggregators\n", "source_name = \"{}/DET/{}CH0:xtdf\" # Source name for raw detector data - filled with karabo_id & module number\n", "\n", "use_dir_creation_date = True # use the creation date of the directory for database time derivation\n", "cal_db_interface = \"tcp://max-exfl016:8015#8025\" # the database interface to use\n", "cal_db_timeout = 300000 # timeout on caldb requests\"\n", "local_output = True # output constants locally\n", "db_output = False # output constants to database\n", "\n", "capacitor_setting = 5 # capacitor_setting for which data was taken \n", "mem_cells = 512 # number of memory cells used\n", "bias_voltage = 250 # detector bias voltage\n", "thresholds_offset_sigma = 3. # bad pixel relative threshold in terms of n sigma offset\n", "thresholds_offset_hard = [400, 1500] # bad pixel hard threshold\n", "thresholds_noise_sigma = 7. # bad pixel relative threshold in terms of n sigma noise\n", "thresholds_noise_hard = [1, 35] # bad pixel hard threshold\n", "\n", "ntrains = 500 # maximum number of trains to use in each gain stage\n", "skip_first_ntrains = 10 # Number of first trains to skip\n", "min_trains = 370 # minimum number of trains needed for each gain stage\n", "\n", "# Parameters for plotting\n", "skip_plots = False # exit after writing corrected files\n", "\n", "high_res_badpix_3d = False # plot bad-pixel summary in high resolution\n", "test_for_normality = False # permorm normality test\n", "inject_cell_order = 'auto' # Include memory cell order as part of the detector condition: auto/always/never\n", "operation_mode = '' # Detector operation mode, optional" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import copy\n", "import multiprocessing\n", "import warnings\n", "from datetime import datetime\n", "from functools import partial\n", "from logging import warning\n", "from pathlib import Path\n", "\n", "warnings.filterwarnings('ignore')\n", "\n", "import dateutil.parser\n", "import matplotlib\n", "import pasha as psh\n", "import scipy.stats\n", "from IPython.display import Latex, Markdown, display\n", "\n", "matplotlib.use(\"agg\")\n", "import matplotlib.patches as patches\n", "import matplotlib.pyplot as plt\n", "\n", "%matplotlib inline\n", "import numpy as np\n", "import tabulate\n", "import yaml\n", "from iCalibrationDB import Conditions, Constants, Detectors, Versions\n", "from XFELDetAna.plotting.heatmap import heatmapPlot\n", "from XFELDetAna.plotting.simpleplot import simplePlot\n", "from extra_data import RunDirectory\n", "\n", "from cal_tools.enums import BadPixels\n", "from cal_tools.lpdlib import make_cell_order_condition\n", "from cal_tools.plotting import (\n", " create_constant_overview,\n", " plot_badpix_3d,\n", " show_overview,\n", " show_processed_modules,\n", ")\n", "from cal_tools.tools import (\n", " get_dir_creation_date,\n", " get_from_db,\n", " get_notebook_name,\n", " get_pdu_from_db,\n", " get_random_db_interface,\n", " get_report,\n", " map_gain_stages,\n", " module_index_to_qm,\n", " parse_runs,\n", " reorder_axes,\n", " run_prop_seq_from_path,\n", " save_const_to_h5,\n", " send_to_db,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gains = np.arange(3)\n", "max_cells = mem_cells\n", "cells = np.arange(max_cells)\n", "gain_names = ['High', 'Medium', 'Low']\n", " \n", "if karabo_da[0] == '-1':\n", " if modules[0] == -1:\n", " modules = list(range(16))\n", " karabo_da = ['LPD{:02d}'.format(i) for i in modules]\n", "else:\n", " modules = [int(x[-2:]) for x in karabo_da]\n", "\n", "capacitor_setting_s = f'{capacitor_setting}pf'\n", "\n", "creation_time = None\n", "if use_dir_creation_date:\n", " creation_time = get_dir_creation_date(in_folder, run_high)\n", " print(\"Using {} as creation time\".format(creation_time))\n", "\n", "if inject_cell_order not in {'auto', 'always', 'never'}:\n", " raise ValueError(\"inject_cell_order must be auto/always/never\")\n", "\n", "run, prop, seq = run_prop_seq_from_path(in_folder)\n", "\n", "cal_db_interface = get_random_db_interface(cal_db_interface)\n", "\n", "display(Markdown('## Evaluated parameters'))\n", "print('CalDB Interface {}'.format(cal_db_interface))\n", "print(\"Proposal: {}\".format(prop))\n", "print(\"Memory cells: {}/{}\".format(mem_cells, max_cells))\n", "print(\"Runs: {}, {}, {}\".format(run_high, run_med, run_low))\n", "print(\"Using DB: {}\".format(db_output))\n", "print(\"Input: {}\".format(in_folder))\n", "print(\"Output: {}\".format(out_folder))\n", "print(\"Bias voltage: {}V\".format(bias_voltage))\n", "print(f\"Capacitor setting: {capacitor_setting_s}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data processing" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "parallel_num_procs = min(6, len(modules)*3)\n", "parallel_num_threads = multiprocessing.cpu_count() // parallel_num_procs\n", "\n", "# the actual characterization\n", "def characterize_module(run_path, channel, gg):\n", " run = RunDirectory(run_path)\n", " det_source = source_name.format(karabo_id, channel)\n", " data = run[det_source, 'image.data'].drop_empty_trains()\n", " data = data[skip_first_ntrains : skip_first_ntrains + ntrains]\n", " cell_ids = run[det_source, 'image.cellId'].drop_empty_trains()\n", " cell_ids = cell_ids[skip_first_ntrains : skip_first_ntrains + ntrains]\n", " \n", " if len(data.train_ids) < min_trains:\n", " raise Exception(f\"Run {run_path} only contains {len(data.train_ids)} trains, but {min_trains} required\")\n", "\n", " im = data.ndarray()\n", " if im.ndim > 3:\n", " im = im[:, 0] # Drop extra dimension\n", " \n", " cellid = cell_ids.ndarray()\n", " cellid_pattern = cell_ids[0].ndarray()\n", " if cellid.ndim > 1:\n", " cellid = cellid[:, 0]\n", " cellid_pattern = cellid_pattern[:, 0]\n", "\n", " # Mask off gain bits, leaving only data\n", " im &= 0b0000111111111111\n", "\n", " im = im.astype(np.float32)\n", " im = reorder_axes(im,\n", " from_order=('frames', 'slow_scan', 'fast_scan'),\n", " to_order=('fast_scan', 'slow_scan', 'frames'),\n", " )\n", "\n", " context = psh.context.ThreadContext(num_workers=parallel_num_threads)\n", " offset = context.alloc(shape=(im.shape[0], im.shape[1], max_cells), dtype=np.float64)\n", " noise = context.alloc(like=offset)\n", " normal_test = context.alloc(like=offset)\n", " def process_cell(worker_id, array_index, cc):\n", " idx = cellid == cc\n", " im_slice = im[..., idx]\n", " if np.any(idx):\n", " offset[..., cc] = np.median(im_slice, axis=2)\n", " noise[..., cc] = np.std(im_slice, axis=2)\n", " if test_for_normality:\n", " _, normal_test[..., cc] = scipy.stats.normaltest(\n", " im[:, :, idx], axis=2)\n", " context.map(process_cell, np.unique(cellid))\n", "\n", " # bad pixels\n", " bp = np.zeros(offset.shape, np.uint32)\n", " # offset related bad pixels\n", " offset_mn = np.nanmedian(offset, axis=(0, 1))\n", " offset_std = np.nanstd(offset, axis=(0, 1))\n", "\n", " bp[(offset < offset_mn-thresholds_offset_sigma*offset_std) |\n", " (offset > offset_mn+thresholds_offset_sigma*offset_std)] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value\n", " bp[(offset < thresholds_offset_hard[0]) | (\n", " offset > thresholds_offset_hard[1])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value\n", " bp[~np.isfinite(offset)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value\n", "\n", " # noise related bad pixels\n", " noise_mn = np.nanmedian(noise, axis=(0, 1))\n", " noise_std = np.nanstd(noise, axis=(0, 1))\n", "\n", " bp[(noise < noise_mn-thresholds_noise_sigma*noise_std) |\n", " (noise > noise_mn+thresholds_noise_sigma*noise_std)] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value\n", " bp[(noise < thresholds_noise_hard[0]) | (\n", " noise > thresholds_noise_hard[1])] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value\n", " bp[~np.isfinite(noise)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value\n", "\n", " idx = (cellid == cellid[0])\n", " return offset, noise, channel, gg, bp, im[12, 12, idx], normal_test, cellid_pattern" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "offset_g = {}\n", "noise_g = {}\n", "badpix_g = {}\n", "data_g = {}\n", "ntest_g = {}\n", "# Should be the same cell order for all modules & all gain stages\n", "cellid_patterns_g = {}\n", "\n", "\n", "inp = []\n", "for gg, run_num in enumerate([run_high, run_med, run_low]):\n", " run_path = Path(in_folder, f\"r{run_num:04d}\")\n", " for i in modules:\n", " inp.append((run_path, i, gg))\n", "\n", "with multiprocessing.Pool(processes=parallel_num_procs) as pool:\n", " results = pool.starmap(characterize_module, inp)\n", "\n", "for ir, r in enumerate(results):\n", " offset, noise, i, gg, bp, data, normal, cellid_pattern = r\n", " if data is None:\n", " warning(f\"No data for module {i} of gain {gg}\")\n", " skip_plots = True\n", " continue\n", " qm = module_index_to_qm(i)\n", " if qm not in offset_g:\n", " offset_g[qm] = np.zeros(offset.shape[:3] + (3,))\n", " print(\"Constant shape:\", offset_g[qm].shape)\n", " noise_g[qm] = np.zeros_like(offset_g[qm])\n", " badpix_g[qm] = np.zeros_like(offset_g[qm], dtype=np.uint32)\n", " data_g[qm] = np.full((ntrains, 3), np.nan)\n", " ntest_g[qm] = np.zeros_like(offset_g[qm])\n", " cellid_patterns_g[qm] = cellid_pattern\n", " else:\n", " if not np.array_equal(cellid_pattern, cellid_patterns_g[qm]):\n", " raise ValueError(\"Inconsistent cell ID pattern between gain stages\")\n", " \n", "\n", " offset_g[qm][..., gg] = offset\n", " noise_g[qm][..., gg] = noise\n", " badpix_g[qm][..., gg] = bp\n", " data_g[qm][:data.shape[0], gg] = data\n", " ntest_g[qm][..., gg] = normal\n", "\n", " hn, cn = np.histogram(data, bins=20)\n", " print(f\"{gain_names[gg]} gain, Module: {qm}. \"\n", " f\"Number of processed trains per cell: {data.shape[0]}.\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if skip_plots:\n", " import sys\n", " print('Skipping plots')\n", " sys.exit(0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Read report path and create file location tuple to add with the injection\n", "proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]\n", "file_loc = 'proposal:{} runs:{} {} {}'.format(proposal, run_low, run_med, run_high)\n", "\n", "report = get_report(metadata_folder)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# TODO: add db_module when received from myMDC\n", "# Create the modules dict of karabo_das and PDUs\n", "qm_dict = {}\n", "for i, k_da in zip(modules, karabo_da):\n", " qm = module_index_to_qm(i)\n", " qm_dict[qm] = {\"karabo_da\": k_da,\n", " \"db_module\": \"\"}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Retrieve existing constants for comparison\n", "clist = [\"Offset\", \"Noise\", \"BadPixelsDark\"]\n", "\n", "dinstance = \"LPD1M1\"\n", "detinst = getattr(Detectors, dinstance)\n", "print('Retrieve pre-existing constants for comparison.')\n", "\n", "old_const = {}\n", "old_mdata = {}\n", "for qm in offset_g.keys():\n", " old_const[qm] = {}\n", " old_mdata[qm] = {}\n", " qm_db = qm_dict[qm]\n", " karabo_da = qm_db[\"karabo_da\"]\n", " cellid_pattern = cellid_patterns_g[qm]\n", "\n", " condition = Conditions.Dark.LPD(\n", " memory_cells=max_cells,\n", " bias_voltage=bias_voltage,\n", " capacitor=capacitor_setting_s,\n", " memory_cell_order=make_cell_order_condition(inject_cell_order, cellid_pattern),\n", " )\n", " for const in clist:\n", " constant = getattr(Constants.LPD, const)()\n", " if not qm_db[\"db_module\"]:\n", " # This should be used in case of running notebook\n", " # by a different method other than myMDC which already\n", " # sends CalCat info.\n", " qm_db[\"db_module\"] = get_pdu_from_db(karabo_id, [karabo_da], constant,\n", " condition, cal_db_interface,\n", " snapshot_at=creation_time)[0]\n", "\n", " data, mdata = get_from_db(karabo_id, karabo_da,\n", " constant,\n", " condition, None,\n", " cal_db_interface,\n", " creation_time=creation_time,\n", " verbosity=2, timeout=cal_db_timeout)\n", "\n", " old_const[qm][const] = data\n", "\n", " if mdata is None or data is None:\n", " old_mdata[qm][const] = {\n", " \"timestamp\": \"Not found\",\n", " \"filepath\": None,\n", " \"h5path\": None\n", " }\n", " else:\n", " timestamp = mdata.calibration_constant_version.begin_at.isoformat()\n", " filepath = Path(\n", " mdata.calibration_constant_version.hdf5path,\n", " mdata.calibration_constant_version.filename\n", " )\n", " h5path = mdata.calibration_constant_version.h5path\n", " old_mdata[qm][const] = {\n", " \"timestamp\": timestamp,\n", " \"filepath\": str(filepath),\n", " \"h5path\": h5path\n", " }\n", "\n", " with open(f\"{out_folder}/module_metadata_{qm}.yml\",\"w\") as fd:\n", " yaml.safe_dump(\n", " {\n", " \"module\": qm,\n", " \"pdu\": qm_db[\"db_module\"],\n", " \"old-constants\": old_mdata[qm]\n", " }, fd)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res = {}\n", "for i in modules:\n", " qm = module_index_to_qm(i)\n", "\n", " res[qm] = {'Offset': offset_g[qm],\n", " 'Noise': noise_g[qm],\n", " 'BadPixelsDark': badpix_g[qm]\n", " }" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Save constants in the calibration DB\n", "md = None\n", "\n", "for qm in res:\n", "\n", " karabo_da = qm_dict[qm][\"karabo_da\"]\n", " db_module = qm_dict[qm][\"db_module\"]\n", " mem_cell_order = make_cell_order_condition(\n", " inject_cell_order, cellid_patterns_g[qm]\n", " )\n", " print(\"Memory cell order:\", mem_cell_order)\n", "\n", " # Do not store empty constants\n", " # In case of 0 trains data_g is initiated with nans and never refilled.\n", " if np.count_nonzero(~np.isnan(data_g[qm]))==0:\n", " print(f\"Constant ({qm}) would be empty, skipping saving\")\n", " continue\n", "\n", " for const in res[qm]:\n", "\n", " dconst = getattr(Constants.LPD, const)()\n", " dconst.data = res[qm][const]\n", "\n", " # set the operating condition\n", "\n", " condition = Conditions.Dark.LPD(memory_cells=max_cells,\n", " bias_voltage=bias_voltage,\n", " capacitor=capacitor_setting_s,\n", " memory_cell_order=mem_cell_order,\n", " )\n", "\n", " if db_output:\n", " md = send_to_db(db_module, karabo_id, dconst, condition,\n", " file_loc, report_path=report,\n", " cal_db_interface=cal_db_interface,\n", " creation_time=creation_time,\n", " timeout=cal_db_timeout)\n", "\n", " if local_output:\n", " md = save_const_to_h5(db_module, karabo_id, dconst, condition,\n", " dconst.data, file_loc, report, creation_time, out_folder)\n", " print(f\"Calibration constant {const} is stored locally.\\n\")\n", "\n", " print(\"Constants parameter conditions are:\\n\")\n", " print(f\"• memory_cells: {max_cells}\\n\"\n", " f\"• bias_voltage: {bias_voltage}\\n\"\n", " f\"• capacitor: {capacitor_setting_s}\\n\"\n", " f\"• memory cell order: {mem_cell_order}\\n\"\n", " f\"• creation_time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\\n\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "show_processed_modules(\n", " dinstance=dinstance,\n", " constants=None,\n", " mnames=[module_index_to_qm(i) for i in modules],\n", " mode=\"position\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Raw pedestal distribution ##\n", "\n", "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", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, grid = plt.subplots(3, 1, sharex=\"col\", sharey=\"row\", figsize=(10, 7))\n", "fig.subplots_adjust(wspace=0, hspace=0)\n", "\n", "for i in modules:\n", " qm = module_index_to_qm(i)\n", " if np.count_nonzero(~np.isnan(data_g[qm])) == 0:\n", " break\n", " for gain in range(3):\n", " data = data_g[qm][:, gain]\n", " offset = np.nanmedian(data)\n", " noise = np.nanstd(data)\n", " xrange = [np.nanmin(data_g[qm]), np.nanmax(data_g[qm])]\n", " if xrange[1] == xrange[0]:\n", " xrange = [0, xrange[0]+xrange[0]//2]\n", " nbins = data_g[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", " grid[gain].hist(data, range=xrange, bins=nbins)\n", " grid[gain].plot([offset-noise, offset-noise], [0, np.nanmax(hn)], \n", " linewidth=1.5, color='red',\n", " label='1 $\\sigma$ deviation')\n", " grid[gain].plot([offset+noise, offset+noise],\n", " [0, np.nanmax(hn)], linewidth=1.5, color='red')\n", " grid[gain].plot([offset, offset], [0, 0],\n", " linewidth=1.5, color='y', label='median')\n", "\n", " grid[gain].plot([np.nanmedian(offset_g[qm][:, :, 12, gain]), \n", " np.nanmedian(offset_g[qm][:, :, 12, gain])],\n", " [0, np.nanmax(hn)], linewidth=1.5, color='green', \n", " label='average over pixels')\n", "\n", " grid[gain].set_xlim(xrange)\n", " grid[gain].set_ylim(0, np.nanmax(hn)*1.1)\n", " grid[gain].set_xlabel(\"Offset value [ADU]\")\n", " grid[gain].set_ylabel(\"# of occurance\")\n", "\n", " if gain == 0:\n", " leg = grid[gain].legend(\n", " loc='upper center', ncol=3, \n", " bbox_to_anchor=(0.1, 0.25, 0.7, 1.0))\n", "\n", " grid[gain].text(820, np.nanmax(hn)*0.4,\n", " \"{} gain\".format(gain_names[gain]), fontsize=20)\n", "\n", " a = plt.axes([.125, .1, 0.775, .8], frame_on=False)\n", " a.patch.set_alpha(0.05)\n", " a.set_xlim(xrange)\n", " plt.plot([offset, offset], [0, 1], linewidth=1.5, color='y')\n", " plt.xticks([])\n", " plt.yticks([])\n", "\n", " ypos = 0.9\n", " x1pos = (np.nanmedian(data_g[qm][:, 0]) +\n", " np.nanmedian(data_g[qm][:, 2]))/2.\n", " x2pos = (np.nanmedian(data_g[qm][:, 2]) +\n", " np.nanmedian(data_g[qm][:, 1]))/2.-10\n", "\n", " plt.annotate(\"\", xy=(np.nanmedian(data_g[qm][:, 0]), ypos), xycoords='data',\n", " xytext=(np.nanmedian(data_g[qm][:, 2]), ypos), textcoords='data',\n", " arrowprops=dict(arrowstyle=\"<->\", connectionstyle=\"arc3\"))\n", "\n", " plt.annotate('{}'.format(np.nanmedian(data_g[qm][:, 0])-np.nanmedian(data_g[qm][:, 2])),\n", " xy=(x1pos, ypos), xycoords='data', xytext=(5, 5), textcoords='offset points')\n", "\n", " plt.annotate(\"\", xy=(np.nanmedian(data_g[qm][:, 2]), ypos), xycoords='data',\n", " xytext=(np.nanmedian(data_g[qm][:, 1]), ypos), textcoords='data',\n", " arrowprops=dict(arrowstyle=\"<->\", connectionstyle=\"arc3\"))\n", "\n", " plt.annotate('{}'.format(np.nanmedian(data_g[qm][:, 2])-np.nanmedian(data_g[qm][:, 1])),\n", " xy=(x2pos, ypos), xycoords='data', xytext=(5, 5), textcoords='offset points')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Normality test ##\n", "\n", "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", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Loop over modules, constants\n", "\n", "if not test_for_normality:\n", " print('Normality test was not requested. Flag `test_for_normality` False')\n", "else:\n", " for i in modules:\n", " qm = module_index_to_qm(i)\n", "\n", " data = np.copy(ntest_g[qm][:,:,:,:])\n", " data[badpix_g[qm][:,:,:,:]>0] = 1.01\n", " \n", " hn,cn = np.histogram(data[:,:,:,0], bins=100)\n", " \n", " d = [{'x': np.arange(100)*0.01+0.01,\n", " 'y': np.histogram(data[:,:,:,0], bins=100)[0],\n", " 'drawstyle': 'steps-pre',\n", " 'label' : 'High gain',\n", " },\n", " {'x': np.arange(100)*0.01+0.01,\n", " 'y': np.histogram(data[:,:,:,1], bins=100)[0],\n", " 'drawstyle': 'steps-pre',\n", " 'label' : 'Medium gain',\n", " },\n", " {'x': np.arange(100)*0.01+0.01,\n", " 'y': np.histogram(data[:,:,:,2], bins=100)[0],\n", " 'drawstyle': 'steps-pre',\n", " 'label' : 'Low gain',\n", " },\n", " ]\n", " \n", "\n", " fig = plt.figure(figsize=(15,15), tight_layout={'pad': 0.5, 'w_pad': 0.3})\n", "\n", " for gain in range(3):\n", " ax = fig.add_subplot(221+gain)\n", " heatmapPlot(data[:,:,12,gain], add_panels=False, cmap='viridis', figsize=(10,10),\n", " y_label='Rows', x_label='Columns',\n", " lut_label='p-Value',\n", " use_axis=ax,\n", " title='p-Value for cell 12, {} gain'.format(gain_names[gain]) )\n", " \n", " ax = fig.add_subplot(224)\n", " _ = simplePlot(d, #aspect=1.6, \n", " x_label = \"p-Value\".format(gain), \n", " y_label=\"# of occurance\",\n", " use_axis=ax,\n", " y_log=False, legend='outside-top-ncol3-frame', legend_pad=0.05, legend_size='5%')\n", " ax.ticklabel_format(style='sci', axis='y', scilimits=(4,6))\n", " " ] }, { "cell_type": "raw", "metadata": {}, "source": [ ".. raw:: latex\n", "\n", " \\newpage" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Single-Cell Overviews ##\n", "\n", "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", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "cell = 12\n", "\n", "for gain in range(3):\n", " display(\n", " Markdown('### Cell-12 overview - {} gain'.format(gain_names[gain])))\n", "\n", " fig = plt.figure(figsize=(18, 22) , tight_layout={'pad': 0.1, 'w_pad': 0.1})\n", " for qm in res:\n", " for iconst, const in enumerate(['Offset', 'Noise', 'BadPixelsDark']):\n", "\n", " ax = fig.add_subplot(321+iconst)\n", "\n", " data = res[qm][const][:, :, 12, gain]\n", " vmax = 1.5 * np.nanmedian(res[qm][const][:, :, 12, gain])\n", " title = const\n", " label = '{} value [ADU]'.format(const)\n", " title = '{} value'.format(const)\n", " if const == 'BadPixelsDark':\n", " vmax = 4\n", " bpix_code = data.astype(np.float32)\n", " bpix_code[bpix_code == 0] = np.nan\n", " title = 'Bad pixel code'\n", " label = title\n", "\n", " cb_labels = ['1 {}'.format(BadPixels.NOISE_OUT_OF_THRESHOLD.name),\n", " '2 {}'.format(BadPixels.OFFSET_NOISE_EVAL_ERROR.name),\n", " '3 {}'.format(BadPixels.OFFSET_OUT_OF_THRESHOLD.name),\n", " '4 {}'.format('MIXED')]\n", "\n", " heatmapPlot(bpix_code, add_panels=False, cmap='viridis',\n", " y_label='Rows', x_label='Columns',\n", " lut_label='', vmax=vmax,\n", " use_axis=ax, cb_ticklabels=cb_labels, cb_ticks = np.arange(4)+1,\n", " title='{}'.format(title))\n", " del bpix_code\n", " else:\n", "\n", " heatmapPlot(data, add_panels=False, cmap='viridis',\n", " y_label='Rows', x_label='Columns',\n", " lut_label=label, vmax=vmax,\n", " use_axis=ax,\n", " title='{}'.format(title))\n", "\n", " for qm in res:\n", " for iconst, const in enumerate(['Offset', 'Noise']):\n", " data = res[qm][const]\n", " dataBP = np.copy(data)\n", " dataBP[res[qm]['BadPixelsDark'] > 0] = -1\n", "\n", " x_ranges = [[0, 1500], [0, 40]]\n", " hn, cn = np.histogram(\n", " data[:, :, :, gain], bins=100, range=x_ranges[iconst])\n", " hnBP, cnBP = np.histogram(dataBP[:, :, :, gain], bins=cn)\n", "\n", " d = [{'x': cn[:-1],\n", " 'y': hn,\n", " 'drawstyle': 'steps-pre',\n", " 'label': 'All data',\n", " },\n", " {'x': cnBP[:-1],\n", " 'y': hnBP,\n", " 'drawstyle': 'steps-pre',\n", " 'label': 'Bad pixels masked',\n", " },\n", " ]\n", "\n", " ax = fig.add_subplot(325+iconst)\n", " _ = simplePlot(d, figsize=(5, 7), aspect=1,\n", " x_label=\"{} value [ADU]\".format(const),\n", " y_label=\"# of occurance\",\n", " title='', legend_pad=0.1, legend_size='10%',\n", " use_axis=ax,\n", " y_log=True, legend='outside-top-2col-frame')\n", "\n", " plt.show()" ] }, { "cell_type": "raw", "metadata": {}, "source": [ ".. raw:: latex\n", "\n", " \\newpage" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cols = {BadPixels.NOISE_OUT_OF_THRESHOLD.value: (BadPixels.NOISE_OUT_OF_THRESHOLD.name, '#FF000080'),\n", " BadPixels.OFFSET_NOISE_EVAL_ERROR.value: (BadPixels.OFFSET_NOISE_EVAL_ERROR.name, '#0000FF80'),\n", " BadPixels.OFFSET_OUT_OF_THRESHOLD.value: (BadPixels.OFFSET_OUT_OF_THRESHOLD.name, '#00FF0080'),\n", " BadPixels.OFFSET_OUT_OF_THRESHOLD.value | BadPixels.NOISE_OUT_OF_THRESHOLD.value: ('MIXED', '#DD00DD80')}\n", "\n", "if high_res_badpix_3d:\n", " display(Markdown(\"\"\"\n", " \n", " ## Global Bad Pixel Behaviour ##\n", "\n", " The following plots shows the results of a bad pixel evaluation for all evaluated memory cells.\n", " Cells are stacked in the Z-dimension, while pixels values in x/y are re-binned with a factor of 2.\n", " This excludes single bad pixels present only in disconnected pixels.\n", " Hence, any bad pixels spanning at least 4 pixels in the x/y-plane, or across at least two memory cells are indicated.\n", " Colors encode the bad pixel type, or mixed type.\n", "\n", " \"\"\"))\n", " # Switch rebin to 1 for full resolution and \n", " # no interpolation for badpixel values.\n", " rebin = 2\n", " for gain in range(3):\n", " display(Markdown('### Bad pixel behaviour - {} gain ###'.format(gain_names[gain])))\n", " for mod, data in badpix_g.items():\n", " plot_badpix_3d(data[...,gain], cols, title='', rebin_fac=rebin)\n", " ax = plt.gca()\n", " leg = ax.get_legend()\n", " leg.set(alpha=0.5)\n", " plt.show()" ] }, { "cell_type": "raw", "metadata": {}, "source": [ ".. raw:: latex\n", "\n", " \\newpage" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary across tiles ##\n", "\n", "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", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "time_summary = []\n", "time_summary.append(f\"The following pre-existing constants are used for comparison:\")\n", "for qm, qm_data in old_mdata.items():\n", " time_summary.append(f\"- Module {qm}\")\n", " for const, const_data in qm_data.items():\n", " time_summary.append(f\" - {const} created at {const_data['timestamp']}\")\n", "display(Markdown(\"\\n\".join(time_summary)))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Loop over modules, constants\n", "for qm in res:\n", " for gain in range(3):\n", " display(Markdown('### Summary across tiles - {} gain'.format(gain_names[gain])))\n", "\n", " for const in res[qm]:\n", " data = np.copy(res[qm][const][:, :, :, gain])\n", "\n", " label = 'Fraction of bad pixels'\n", "\n", " if const != 'BadPixelsDark':\n", " data[badpix_g[qm][:, :, :, gain] > 0] = np.nan\n", " label = '{} value [ADU]'.format(const)\n", " else:\n", " data[data>0] = 1.0\n", "\n", " data = data.reshape(\n", " int(data.shape[0] / 32),\n", " 32,\n", " int(data.shape[1] / 128),\n", " 128,\n", " data.shape[2])\n", " data = np.nanmean(data, axis=(1, 3)).swapaxes(\n", " 0, 2).reshape(512, 16)\n", "\n", " fig = plt.figure(figsize=(15, 6))\n", " ax = fig.add_subplot(121)\n", "\n", " _ = heatmapPlot(data[:510, :], add_panels=True,\n", " y_label='Momery Cell ID', x_label='Tile ID',\n", " lut_label=label, use_axis=ax,\n", " panel_y_label=label, panel_x_label=label,\n", " cmap='viridis', # cb_loc='right',cb_aspect=15,\n", " x_ticklabels=np.arange(16)+1,\n", " x_ticks=np.arange(16)+0.5)\n", "\n", " if old_const[qm][const] is not None:\n", " ax = fig.add_subplot(122)\n", "\n", " dataold = np.copy(old_const[qm][const][:, :, :, gain])\n", "\n", " label = '$\\Delta$ {}'.format(label)\n", "\n", " if const != 'BadPixelsDark':\n", " if old_const[qm]['BadPixelsDark'] is not None:\n", " dataold[old_const[qm]['BadPixelsDark'][:, :, :, gain] > 0] = np.nan\n", " else:\n", " dataold[:] = np.nan\n", " else:\n", " dataold[dataold>0]=1.0\n", "\n", " dataold = dataold.reshape(\n", " int(dataold.shape[0] / 32),\n", " 32,\n", " int(dataold.shape[1] / 128),\n", " 128,\n", " dataold.shape[2])\n", " dataold = np.nanmean(dataold, axis=(\n", " 1, 3)).swapaxes(0, 2).reshape(512, 16)\n", " dataold = dataold - data\n", "\n", " _ = heatmapPlot(dataold[:510, :], add_panels=True,\n", " y_label='Momery Cell ID', x_label='Tile ID',\n", " lut_label=label, use_axis=ax,\n", " panel_y_label=label, panel_x_label=label,\n", " cmap='viridis', # cb_loc='right',cb_aspect=15,\n", " x_ticklabels=np.arange(16)+1,\n", " x_ticks=np.arange(16)+0.5)\n", " plt.show()" ] }, { "cell_type": "raw", "metadata": {}, "source": [ ".. raw:: latex\n", "\n", " \\newpage" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Variation of offset and noise across Tiles and ASICs ##\n", "\n", "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", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Loop over modules, constants\n", "for qm in res:\n", " for gain in range(3):\n", " display(Markdown('### Variation of offset and noise across ASICs - {} gain'.format(gain_names[gain])))\n", "\n", " fig = plt.figure(figsize=(15, 6))\n", " for iconst, const in enumerate(['Offset', 'Noise']):\n", " data = np.copy(res[qm][const][:, :, :, gain])\n", " data[badpix_g[qm][:, :, :, gain] > 0] = np.nan\n", " label = '$\\sigma$ {} [ADU]'.format(const)\n", "\n", " dataA = np.nanmean(data, axis=2) # average over cells\n", " dataA = dataA.reshape(8, 32, 16, 16)\n", " dataA = np.nanstd(dataA, axis=(0, 2)) # average across ASICs\n", "\n", " ax = fig.add_subplot(121+iconst)\n", " _ = heatmapPlot(dataA, add_panels=True,\n", " y_label='rows', x_label='columns',\n", " lut_label=label, use_axis=ax,\n", " panel_y_label=label, panel_x_label=label,\n", " cmap='viridis'\n", " )\n", "\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Loop over modules, constants\n", "for qm in res:\n", " for gain in range(3):\n", " display(Markdown('### Variation of offset and noise across tiles - {} gain'.format(gain_names[gain])))\n", "\n", " fig = plt.figure(figsize=(15, 6))\n", " for iconst, const in enumerate(['Offset', 'Noise']):\n", " data = np.copy(res[qm][const][:, :, :, gain])\n", " data[badpix_g[qm][:, :, :, gain] > 0] = np.nan\n", " label = '$\\sigma$ {} [ADU]'.format(const)\n", "\n", " dataT = data.reshape(\n", " int(data.shape[0] / 32),\n", " 32,\n", " int(data.shape[1] / 128),\n", " 128,\n", " data.shape[2])\n", " dataT = np.nanstd(dataT, axis=(0, 2))\n", " dataT = np.nanmean(dataT, axis=2)\n", "\n", " ax = fig.add_subplot(121+iconst)\n", " _ = heatmapPlot(dataT, add_panels=True,\n", " y_label='rows', x_label='columns',\n", " lut_label=label, use_axis=ax,\n", " panel_y_label=label, panel_x_label=label,\n", " cmap='viridis')\n", " plt.show()" ] }, { "cell_type": "raw", "metadata": {}, "source": [ ".. raw:: latex\n", "\n", " \\newpage" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aggregate values and per cell behaviour ##\n", "\n", "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", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Loop over modules, constants\n", "for qm in res:\n", " for gain in range(3):\n", " display(Markdown('### Mean over pixels - {} gain'.format(gain_names[gain])))\n", "\n", " fig = plt.figure(figsize=(9,11))\n", "\n", " for iconst, const in enumerate(res[qm]):\n", "\n", " ax = fig.add_subplot(311+iconst)\n", "\n", " data = res[qm][const][:,:,:510,gain]\n", " if const == 'BadPixelsDark':\n", " data[data>0] = 1.0\n", "\n", " dataBP = np.copy(data)\n", " dataBP[badpix_g[qm][:,:,:510,gain]>0] = -10\n", "\n", " data = np.nanmean(data, axis=(0,1))\n", " dataBP = np.nanmean(dataBP, axis=(0,1))\n", "\n", " d = [{'y': data,\n", " 'x': np.arange(data.shape[0]),\n", " 'drawstyle': 'steps-mid',\n", " 'label' : 'All data'\n", " }\n", " ]\n", "\n", " if const != 'BadPixelsDark':\n", " d.append({'y': dataBP,\n", " 'x': np.arange(data.shape[0]),\n", " 'drawstyle': 'steps-mid',\n", " 'label' : 'good pixels only'\n", " })\n", " y_title = \"{} value [ADU]\".format(const)\n", " title = \"{} value, {} gain\".format(const, gain_names[gain])\n", " else:\n", " y_title = \"Fraction of Bad Pixels\"\n", " title = \"Fraction of Bad Pixels, {} gain\".format(gain_names[gain])\n", "\n", " data_min = np.min([data, dataBP])if const != 'BadPixelsDark' else np.min([data])\n", " data_max = np.max([data[20:], dataBP[20:]])\n", " data_dif = data_max - data_min\n", "\n", " local_max = np.max([data[200:300], dataBP[200:300]])\n", " frac = 0.35\n", " new_max = (local_max - data_min*(1-frac))/frac\n", " new_max = np.max([data_max, new_max])\n", "\n", " _ = simplePlot(d, figsize=(10,10), aspect=2, xrange=(-12, 510),\n", " x_label = 'Memory Cell ID', \n", " y_label=y_title, use_axis=ax,\n", " title=title,\n", " title_position=[0.5, 1.15], \n", " inset='xy-coord-right', inset_x_range=(0,20), inset_indicated=True,\n", " inset_labeled=True, inset_coord=[0.2,0.5,0.6,0.95],\n", " inset_lw = 1.0, y_range = [data_min-data_dif*0.05, new_max+data_dif*0.05],\n", " y_log=False, legend='outside-top-ncol2-frame', legend_size='18%',\n", " legend_pad=0.00)\n", "\n", " plt.tight_layout(pad=1.08, h_pad=0.35)\n", "\n", " plt.show()" ] }, { "cell_type": "raw", "metadata": {}, "source": [ ".. raw:: latex\n", "\n", " \\newpage" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary tables ##\n", "\n", "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", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "table = []\n", "bits = [\n", " BadPixels.NOISE_OUT_OF_THRESHOLD,\n", " BadPixels.OFFSET_OUT_OF_THRESHOLD,\n", " BadPixels.OFFSET_NOISE_EVAL_ERROR\n", "]\n", "\n", "for qm in res:\n", " for gain in range(3):\n", "\n", " l_data = []\n", " l_data_old = []\n", "\n", " data = np.copy(res[qm]['BadPixelsDark'][:,:,:,gain])\n", " l_data.append(len(data[data>0].flatten()))\n", " for bit in bits:\n", " l_data.append(np.count_nonzero(badpix_g[qm][:,:,:,gain] & bit.value))\n", "\n", " if old_const[qm]['BadPixelsDark'] is not None:\n", " old_const[qm]['BadPixelsDark'] = old_const[qm]['BadPixelsDark'].astype(np.uint32)\n", " dataold = np.copy(old_const[qm]['BadPixelsDark'][:, :, :, gain])\n", " l_data_old.append(len(dataold[dataold>0].flatten()))\n", " for bit in bits:\n", " l_data_old.append(np.count_nonzero(old_const[qm]['BadPixelsDark'][:, :, :, gain] & bit.value))\n", "\n", " l_data_name = ['All bad pixels', 'NOISE_OUT_OF_THRESHOLD', \n", " 'OFFSET_OUT_OF_THRESHOLD', 'OFFSET_NOISE_EVAL_ERROR']\n", "\n", " l_threshold = ['', f'{thresholds_noise_sigma}', f'{thresholds_offset_sigma}',\n", " f'{thresholds_offset_hard}/{thresholds_noise_hard}']\n", "\n", " for i in range(len(l_data)):\n", " line = [f'{l_data_name[i]}, gain {gain_names[gain]}', l_threshold[i], l_data[i]]\n", "\n", " if old_const[qm]['BadPixelsDark'] is not None:\n", " line += [l_data_old[i]]\n", " else:\n", " line += ['-']\n", "\n", " table.append(line)\n", " table.append(['', '', '', ''])\n", "\n", "display(Markdown('''\n", "\n", "### Number of bad pixels ###\n", "\n", "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.\n", "\n", "'''))\n", "if len(table)>0:\n", " md = display(Latex(tabulate.tabulate(table, tablefmt='latex', \n", " headers=[\"Pixel type\", \"Threshold\",\n", " \"New constant\", \"Old constant\"]))) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "header = ['Parameter', \n", " \"New constant\", \"Old constant \", \n", " \"New constant\", \"Old constant \", \n", " \"New constant\", \"Old constant \"]\n", "\n", "for const in ['Offset', 'Noise']:\n", " table = [['','High gain', 'High gain', 'Medium gain', 'Medium gain', 'Low gain', 'Low gain']]\n", " for qm in res:\n", "\n", " data = np.copy(res[qm][const])\n", " data[res[qm]['BadPixelsDark']>0] = np.nan\n", "\n", " if old_const[qm][const] is not None and old_const[qm]['BadPixelsDark'] is not None :\n", " dataold = np.copy(old_const[qm][const])\n", " dataold[old_const[qm]['BadPixelsDark']>0] = np.nan\n", "\n", " f_list = [np.nanmedian, np.nanmean, np.nanstd, np.nanmin, np.nanmax]\n", " n_list = ['Median', 'Mean', 'Std', 'Min', 'Max']\n", "\n", " for i, f in enumerate(f_list):\n", " line = [n_list[i]]\n", " for gain in range(3):\n", " line.append('{:6.1f}'.format(f(data[...,gain])))\n", " if old_const[qm][const] is not None and old_const[qm]['BadPixelsDark'] is not None:\n", " line.append('{:6.1f}'.format(f(dataold[...,gain])))\n", " else:\n", " line.append('-')\n", "\n", " table.append(line)\n", "\n", " display(Markdown('### {} [ADU], good pixels only ###'.format(const)))\n", " md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=header))) " ] } ], "metadata": { "kernelspec": { "display_name": "Offline Cal", "language": "python", "name": "offline-cal" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 2 }