{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# AGIPD Characterize Dark Images #\n", "\n", "Author: European XFEL Detector Group, Version: 2.0\n", "\n", "The following code analyzes a set of dark images taken with the AGIPD detector to deduce detector offsets , noise, bad-pixel maps and thresholding. All four types of constants are evaluated per-pixel and per-memory cell. Data for the detector's three gain stages needs to be present, separated into separate runs.\n", "\n", "The evaluated calibration constants are stored locally and injected in the calibration data base.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "in_folder = \"/gpfs/exfel/d/raw/CALLAB/202031/p900113\" # path to input data, required\n", "out_folder = \"\" # path to output to, required\n", "modules = [-1] # list of modules to evaluate, RANGE ALLOWED\n", "run_high = 9985 # run number in which high gain data was recorded, required\n", "run_med = 9984 # run number in which medium gain data was recorded, required\n", "run_low = 9983 # run number in which low gain data was recorded, required\n", "operation_mode = \"ADAPTIVE_GAIN\" # Detector operation mode, optional (defaults to \"ADAPTIVE_GAIN\")\n", "\n", "karabo_id = \"HED_DET_AGIPD500K2G\" # karabo karabo_id\n", "karabo_da = ['-1'] # a list of data aggregators names, Default [-1] for selecting all data aggregators\n", "receiver_template = \"{}CH0\" # inset for receiver devices\n", "instrument_source_template = '{}/DET/{}:xtdf' # path in the HDF5 file to images\n", "ctrl_source_template = '{}/MDL/FPGA_COMP' # path to control information\n", "karabo_id_control = \"HED_EXP_AGIPD500K2G\" # karabo-id for control device '\n", "\n", "use_dir_creation_date = True # use dir creation date as data production reference date\n", "cal_db_interface = \"tcp://max-exfl016:8020\" # the database interface to use\n", "cal_db_timeout = 3000000 # timeout on caldb requests\"\n", "local_output = True # output constants locally\n", "db_output = False # output constants to database\n", "\n", "mem_cells = 0 # number of memory cells used, set to 0 to automatically infer\n", "bias_voltage = 0 # bias voltage, set to 0 to use stored value in slow data.\n", "gain_setting = -1 # the gain setting, use -1 to use value stored in slow data.\n", "gain_mode = -1 # gain mode, use -1 to use value stored in slow data.\n", "integration_time = -1 # integration time, negative values for auto-detection.\n", "acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine\n", "interlaced = False # assume interlaced data format, for data prior to Dec. 2017\n", "\n", "thresholds_offset_sigma = 3. # offset sigma thresholds for offset deduced bad pixels\n", "thresholds_offset_hard = [0, 0] # For setting the same threshold offset for the 3 gains. Left for backcompatability. Default [0, 0] to take the following parameters.\n", "thresholds_offset_hard_hg = [3000, 7000] # High-gain thresholds in absolute ADU terms for offset deduced bad pixels\n", "thresholds_offset_hard_mg = [6000, 10000] # Medium-gain thresholds in absolute ADU terms for offset deduced bad pixels\n", "thresholds_offset_hard_lg = [6000, 10000] # Low-gain thresholds in absolute ADU terms for offset deduced bad pixels\n", "thresholds_offset_hard_hg_fixed = [3500, 6500] # Same as thresholds_offset_hard_hg, but for fixed gain operation\n", "thresholds_offset_hard_mg_fixed = [3500, 6500] # Same as thresholds_offset_hard_mg, but for fixed gain operation\n", "thresholds_offset_hard_lg_fixed = [3500, 6500] # Same as thresholds_offset_hard_lg, but for fixed gain operation\n", "\n", "thresholds_noise_sigma = 5. # noise sigma thresholds for offset deduced bad pixels\n", "thresholds_noise_hard = [0, 0] # For setting the same threshold noise for the 3 gains. Left for backcompatability. Default [0, 0] to take the following parameters.\n", "thresholds_noise_hard_hg = [4, 20] # High-gain thresholds in absolute ADU terms for offset deduced bad pixels\n", "thresholds_noise_hard_mg = [4, 20] # Medium-gain thresholds in absolute ADU terms for offset deduced bad pixels\n", "thresholds_noise_hard_lg = [4, 20] # Low-gain thresholds in absolute ADU terms for offset deduced bad pixels\n", "\n", "thresholds_gain_sigma = 5. # Gain separation sigma threshold\n", "max_trains = 0 # Maximum number of trains to use for processing dark. Set to 0 to process all available trains.\n", "min_trains = 1 # Miniumum number of trains for processing dark. If raw folder has less than minimum trains processing is stopped.\n", "high_res_badpix_3d = False # set this to True if you need high-resolution 3d bad pixel plots. ~7mins extra time for 64 memory cells\n", "\n", "# This is used if modules is not specified:\n", "def find_modules(in_folder, run_high, modules):\n", " if (modules is not None) and modules != [-1]:\n", " return modules\n", " from pathlib import Path\n", " import re\n", " modules = set()\n", " for file in Path(in_folder, f'r{run_high:04d}').iterdir():\n", " m = re.search(r'-AGIPD(\\d{2})-', file.name)\n", " if m:\n", " modules.add(int(m.group(1)))\n", " return sorted(modules)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import itertools\n", "import multiprocessing\n", "import os\n", "from collections import OrderedDict\n", "from datetime import timedelta\n", "from pathlib import Path\n", "from typing import List, Tuple\n", "\n", "import dateutil.parser\n", "import matplotlib\n", "import numpy as np\n", "import pasha as psh\n", "import psutil\n", "import tabulate\n", "import yaml\n", "from IPython.display import Latex, Markdown, display\n", "from extra_data import RunDirectory\n", "\n", "matplotlib.use('agg')\n", "\n", "import iCalibrationDB\n", "import matplotlib.pyplot as plt\n", "from cal_tools.agipdlib import AgipdCtrl\n", "from cal_tools.enums import AgipdGainMode, BadPixels\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_pdu_from_db,\n", " get_random_db_interface,\n", " get_report,\n", " map_gain_stages,\n", " module_index_to_qm,\n", " run_prop_seq_from_path,\n", " save_const_to_h5,\n", " send_to_db,\n", ")\n", "\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# insert control device if format string (does nothing otherwise)\n", "ctrl_src = ctrl_source_template.format(karabo_id_control)\n", "\n", "runs_dict = OrderedDict()\n", "\n", "for gain_idx, (run_name, run_number) in enumerate(zip(\n", " [\"high\", \"med\", \"low\"],\n", " [run_high, run_med, run_low]\n", ")):\n", " runs_dict[run_name] = {\n", " \"number\": run_number,\n", " \"gain\": gain_idx,\n", " \"dc\": RunDirectory(f'{in_folder}/r{run_number:04d}/')\n", " }\n", "\n", "creation_time=None\n", "if use_dir_creation_date:\n", " creation_time = get_dir_creation_date(in_folder, run_high)\n", "\n", "print(f\"Using {creation_time} as creation time of constant.\")\n", "\n", "run, prop, seq = run_prop_seq_from_path(in_folder)\n", "\n", "# Read report path and create file location tuple to add with the injection\n", "file_loc = f\"proposal:{prop} runs:{run_low} {run_med} {run_high}\"\n", "\n", "report = get_report(out_folder)\n", "cal_db_interface = get_random_db_interface(cal_db_interface)\n", "print(f'Calibration database interface: {cal_db_interface}')\n", "\n", "instrument = karabo_id.split(\"_\")[0]\n", "\n", "if instrument == \"SPB\":\n", " dinstance = \"AGIPD1M1\"\n", " nmods = 16\n", "elif instrument == \"MID\":\n", " dinstance = \"AGIPD1M2\"\n", " nmods = 16\n", "elif instrument == \"HED\":\n", " dinstance = \"AGIPD500K\"\n", " nmods = 8\n", "\n", "instrument_src = instrument_source_template.format(karabo_id, receiver_template)\n", "run_numbers = [run_high, run_med, run_low]\n", "\n", "def create_karabo_da_list(modules):\n", " return([\"AGIPD{:02d}\".format(i) for i in modules])\n", "\n", "if karabo_da[0] == '-1':\n", " if modules[0] == -1:\n", " modules = list(range(nmods))\n", " karabo_da = create_karabo_da_list(modules)\n", "else:\n", " modules = [int(x[-2:]) for x in karabo_da]\n", "\n", "print(f\"Detector in use is {karabo_id}\")\n", "print(f\"Instrument {instrument}\")\n", "print(f\"Detector instance {dinstance}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create out_folder if it doesn't exist.\n", "Path(out_folder).mkdir(parents=True, exist_ok=True)\n", "\n", "n_files = 0\n", "total_file_sizes = 0\n", "max_trains_list = []\n", "\n", "for run_dict in runs_dict.values():\n", " missing_modules = []\n", " image_dc = run_dict[\"dc\"].select(f\"{karabo_id_control}*\", \"*\", require_all=True)\n", " # This is important in case of no slurm parallelization over modules is done.\n", " # (e.g. running notebook interactively)\n", " sources_l = [(f\"{karabo_id_control}*\", \"*\")]\n", " sources_l += [(instrument_src.format(m), \"*\") for m in modules]\n", " image_dc = run_dict[\"dc\"].select(sources_l, require_all=True)\n", " # validate that there are trains and that data sources are\n", " # present for any of the selected modules.\n", " if (\n", " len(image_dc.train_ids) == 0 or\n", " not np.any([\n", " karabo_id in s for s in run_dict[\"dc\"].select(sources_l, require_all=True).all_sources]) # noqa\n", " ):\n", " raise ValueError(f\"No images to process for run: {run_dict['number']}\")\n", "\n", " max_trains_list.append(len(image_dc.train_ids))\n", "\n", " # update run_dc with selected module sources\n", " run_dict[\"dc\"] = image_dc\n", "\n", "# Update modules and karabo_da lists based on available modules to processes.\n", "modules = [m for m in modules if m not in missing_modules]\n", "karabo_da = create_karabo_da_list(modules)\n", "\n", "# Remodifing run data collections to display actual total files number and size. \n", "for run_dict in runs_dict.values():\n", " file_sizes = [os.path.getsize(f.filename) / 1e9 for f in run_dict[\"dc\"].deselect(f\"{karabo_id_control}*\").files]\n", " total_file_sizes += sum(file_sizes)\n", " n_files += len(file_sizes)\n", "\n", "print(f\"Will process data in a total of {n_files} files ({total_file_sizes:.02f} GB).\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read and validate the runs control data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def read_run_conditions(runs_dict: dict):\n", " agipd_cond = AgipdCtrl(\n", " run_dc=runs_dict[\"dc\"],\n", " image_src=instrument_src_mod,\n", " ctrl_src=ctrl_src,\n", " )\n", " cond_dict[\"runs\"].append(runs_dict[\"number\"])\n", " if acq_rate == 0:\n", " cond_dict[\"acq_rate\"].append(agipd_cond.get_acq_rate())\n", " if mem_cells == 0:\n", " cond_dict[\"mem_cells\"].append(agipd_cond.get_num_cells())\n", " if gain_setting == -1: \n", " cond_dict[\"gain_setting\"].append(\n", " agipd_cond.get_gain_setting(creation_time))\n", " if bias_voltage == 0.:\n", " cond_dict[\"bias_voltage\"].append(\n", " agipd_cond.get_bias_voltage(karabo_id_control))\n", " if integration_time == -1:\n", " cond_dict[\"integration_time\"].append(\n", " agipd_cond.get_integration_time())\n", " if gain_mode == -1:\n", " cond_dict[\"gain_mode\"].append(agipd_cond.get_gain_mode())\n", " else:\n", " cond_dict[\"gain_mode\"].append(AgipdGainMode(gain_mode))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def validate_gain_modes(gain_modes: List[AgipdGainMode]):\n", " # Validate that gain modes are not a mix of adaptive and fixed gain.\n", " if all(\n", " gm == AgipdGainMode.ADAPTIVE_GAIN for gm in gain_modes\n", " ):\n", " fixed_gain_mode = False\n", " elif any(\n", " gm == AgipdGainMode.ADAPTIVE_GAIN for gm in gain_modes\n", " ):\n", " raise ValueError(\n", " f\"ERROR: Given runs {self.read_conditions['run_number']}\"\n", " \" have a mix of ADAPTIVE and FIXED gain modes: \"\n", " f\"{self.read_conditions['gain_mode']}.\"\n", " )\n", " else:\n", " fixed_gain_mode = True\n", " return fixed_gain_mode" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Read slow data from 1st channel only.\n", "# Read all modules in one notebook and validate the conditions across detectors?\n", "# Currently slurm jobs run per one module.\n", "\n", "# TODO: what if first module is not available. Maybe only channel 2 available\n", "instrument_src_mod = instrument_src.format(modules[0])\n", "\n", "cond_dict = dict()\n", "fixed_gain_mode = None\n", "\n", "with multiprocessing.Manager() as manager:\n", " cond_dict[\"runs\"] = manager.list()\n", " cond_dict[\"acq_rate\"] = manager.list()\n", " cond_dict[\"mem_cells\"] = manager.list()\n", " cond_dict[\"gain_setting\"] = manager.list()\n", " cond_dict[\"gain_mode\"] = manager.list()\n", " cond_dict[\"bias_voltage\"] = manager.list()\n", " cond_dict[\"integration_time\"] = manager.list()\n", "\n", " with multiprocessing.Pool(processes=len(modules)) as pool:\n", " pool.starmap(read_run_conditions, zip(runs_dict.values()))\n", "\n", " for cond, vlist in cond_dict.items():\n", " if cond == \"runs\":\n", " continue\n", " elif cond == \"gain_mode\":\n", " fixed_gain_mode = validate_gain_modes(cond_dict[\"gain_mode\"])\n", " if not all(x == vlist[0] for x in vlist):\n", " # TODO: raise ERROR??\n", " print(\n", " f\"WARNING: {cond} is not the same for the runs \"\n", " f\"{cond_dict['runs']} with values\"\n", " f\" of {cond_dict[cond]}, respectively.\"\n", " )\n", " if cond_dict[\"acq_rate\"]: acq_rate = cond_dict[\"acq_rate\"][0]\n", " if cond_dict[\"mem_cells\"]: mem_cells = cond_dict[\"mem_cells\"][0]\n", " if cond_dict[\"gain_setting\"]: gain_setting = cond_dict[\"gain_setting\"][0]\n", " if cond_dict[\"gain_mode\"]: gain_mode = list(cond_dict[\"gain_mode\"])\n", " if cond_dict[\"bias_voltage\"]: bias_voltage = cond_dict[\"bias_voltage\"][0]\n", " if cond_dict[\"integration_time\"]: integration_time = cond_dict[\"integration_time\"][0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Determine the gain operation mode based on the gain_mode stored in control h5file.\n", "if operation_mode not in (\"ADAPTIVE_GAIN\", \"FIXED_GAIN\"):\n", " print(f\"WARNING: unknown operation_mode \\\"{operation_mode}\\\" parameter set\")\n", "\n", "if (\n", " gain_mode == [\n", " AgipdGainMode.FIXED_HIGH_GAIN,\n", " AgipdGainMode.FIXED_MEDIUM_GAIN,\n", " AgipdGainMode.FIXED_LOW_GAIN\n", " ] and\n", " operation_mode == \"ADAPTIVE_GAIN\"\n", "):\n", " print(\n", " \"WARNING: operation_mode parameter is ADAPTIVE_GAIN, \"\n", " \"slow data indicates FIXED_GAIN.\")\n", "elif not fixed_gain_mode and operation_mode == \"FIXED_GAIN\":\n", " print(\n", " \"WARNING: operation_mode parameter is FIXED_GAIN, \"\n", " \"slow data indicates ADAPTIVE_GAIN\")\n", "elif not all(gm == AgipdGainMode.ADAPTIVE_GAIN for gm in gain_mode):\n", " raise ValueError(\n", " \"ERROR: Wrong arrangment of given dark runs. \"\n", " f\"Given runs' gain_modes are {gain_mode} for runs: {runs}.\"\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Parameters are:\")\n", "print(f\"Proposal: {prop}\")\n", "print(f\"Acquisition rate: {acq_rate}\")\n", "print(f\"Memory cells: {mem_cells}\")\n", "print(f\"Runs: {run_numbers}\")\n", "print(f\"Interlaced mode: {interlaced}\")\n", "print(f\"Using DB: {db_output}\")\n", "print(f\"Input: {in_folder}\")\n", "print(f\"Output: {out_folder}\")\n", "print(f\"Bias voltage: {bias_voltage}V\")\n", "print(f\"Gain setting: {gain_setting}\")\n", "print(f\"Integration time: {integration_time}\")\n", "print(f\"Operation mode is {'fixed' if fixed_gain_mode else 'adaptive'} gain mode\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if thresholds_offset_hard != [0, 0]:\n", " # if set, this will override the individual parameters\n", " thresholds_offset_hard = [thresholds_offset_hard] * 3\n", "elif fixed_gain_mode:\n", " thresholds_offset_hard = [\n", " thresholds_offset_hard_hg_fixed,\n", " thresholds_offset_hard_mg_fixed,\n", " thresholds_offset_hard_lg_fixed,\n", " ]\n", "else:\n", " thresholds_offset_hard = [\n", " thresholds_offset_hard_hg,\n", " thresholds_offset_hard_mg,\n", " thresholds_offset_hard_lg,\n", " ]\n", "print(\"Will use the following hard offset thresholds\")\n", "for name, value in zip((\"High\", \"Medium\", \"Low\"), thresholds_offset_hard):\n", " print(f\"- {name} gain: {value}\")\n", "\n", "if thresholds_noise_hard != [0, 0]:\n", " thresholds_noise_hard = [thresholds_noise_hard] * 3\n", "else:\n", " thresholds_noise_hard = [\n", " thresholds_noise_hard_hg,\n", " thresholds_noise_hard_mg,\n", " thresholds_noise_hard_lg,\n", " ]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Check if max_trains can be processed.\n", "\n", "# more relevant if running on multiple modules (i.e. within notebook)\n", "# mem_cells * gains * n_constants * modules * agipd_[x,y]image_size * 2\n", "av_mem = psutil.virtual_memory().available\n", "possible_trains = av_mem // (352 * 3 * 3 * len(modules) * 131072 * 2)\n", "if max_trains == 0:\n", " max_trains = max(max_trains_list)\n", "if max_trains > possible_trains:\n", " max_trains = possible_trains\n", " print(\n", " f\"WARNING: available memory for processing is { av_mem / 1e9:.02f} GB.\"\n", " f\" Modifing max_trains to process to {max_trains}\")\n", "\n", "for run_dict in runs_dict.values():\n", " run_dict[\"dc\"] = run_dict[\"dc\"].select_trains(np.s_[:max_trains])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate Offsets, Noise and Thresholds ##\n", "\n", "The calculation is performed per-pixel and per-memory-cell. Offsets are simply the median value for a set of dark data taken at a given gain, noise the standard deviation, and gain-bit values the medians of the gain array." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "parallel_num_procs = min(12, len(modules)*3)\n", "parallel_num_threads = multiprocessing.cpu_count() // parallel_num_procs\n", "print(f\"Will use {parallel_num_procs} processes with {parallel_num_threads} threads each\")\n", "\n", "def characterize_module(\n", " channel: int, runs_dict: dict,\n", ") -> Tuple[int, int, np.array, np.array, np.array, np.array, np.array]:\n", "\n", " # Select the corresponding module channel.\n", " instrument_src_mod = instrument_src.format(channel)\n", "\n", " run_dc = runs_dict[\"dc\"]\n", " gain_index = runs_dict[\"gain\"]\n", "\n", " if run_dc[instrument_src_mod, \"image.data\"].shape[0] < min_trains:\n", " print(\n", " f\"WARNING: {run_dc.files} have less than \"\n", " \"minimum trains: {min_trains}.\")\n", "\n", " # Read module's image and cellId data.\n", " im = run_dc[instrument_src_mod, \"image.data\"].ndarray()\n", " cell_ids = np.squeeze(run_dc[instrument_src_mod, \"image.cellId\"].ndarray())\n", "\n", " local_thresholds_offset_hard = thresholds_offset_hard[gain_index]\n", " local_thresholds_noise_hard = thresholds_noise_hard[gain_index] \n", "\n", " if interlaced:\n", " if not fixed_gain_mode:\n", " ga = im[1::2, 0, ...]\n", " im = im[0::2, 0, ...].astype(np.float32)\n", " cell_ids = cell_ids[::2]\n", " else:\n", " if not fixed_gain_mode:\n", " ga = im[:, 1, ...]\n", " im = im[:, 0, ...].astype(np.float32)\n", " im = np.transpose(im)\n", " if not fixed_gain_mode:\n", " ga = np.transpose(ga)\n", "\n", " context = psh.context.ThreadContext(num_workers=parallel_num_threads)\n", " offset = context.alloc(shape=(im.shape[0], im.shape[1], mem_cells), dtype=np.float64)\n", " noise = context.alloc(like=offset)\n", "\n", " if fixed_gain_mode:\n", " gains = None\n", " gains_std = None\n", " else:\n", " gains = context.alloc(like=offset)\n", " gains_std = context.alloc(like=offset)\n", "\n", " def process_cell(worker_id, array_index, cell_number):\n", " cell_slice_index = (cell_ids == cell_number)\n", " im_slice = im[..., cell_slice_index]\n", " offset[..., cell_number] = np.median(im_slice, axis=2)\n", " noise[..., cell_number] = np.std(im_slice, axis=2)\n", " if not fixed_gain_mode:\n", " ga_slice = ga[..., cell_slice_index]\n", " gains[..., cell_number] = np.median(ga_slice, axis=2)\n", " gains_std[..., cell_number] = np.std(ga_slice, axis=2)\n", " context.map(process_cell, np.unique(cell_ids))\n", "\n", " # bad pixels\n", " bp = np.zeros_like(offset, dtype=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\n", " bp[(offset < local_thresholds_offset_hard[0]) |\n", " (offset > local_thresholds_offset_hard[1])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD\n", " bp[~np.isfinite(offset)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR\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", " bp[(noise < noise_mn-thresholds_noise_sigma*noise_std) |\n", " (noise > noise_mn+thresholds_noise_sigma*noise_std)] |= BadPixels.NOISE_OUT_OF_THRESHOLD\n", " bp[(noise < local_thresholds_noise_hard[0]) | (noise > local_thresholds_noise_hard[1])] |= BadPixels.NOISE_OUT_OF_THRESHOLD\n", " bp[~np.isfinite(noise)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR\n", "\n", " return channel, gain_index, offset, noise, gains, gains_std, bp" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "with multiprocessing.Pool(processes=parallel_num_procs) as pool:\n", " results = pool.starmap(\n", " characterize_module, itertools.product(modules, list(runs_dict.values())))\n", "\n", "# mapped values for processing 2 modules example:\n", "# [\n", "# 0, {\"gain\": 0, \"run_number\": <run-high>, \"dc\": <high-dc>},\n", "# 0, {\"gain\": 1, \"run_number\": <run-med>, \"dc\": <med-dc>},\n", "# 0, {\"gain\": 2, \"run_number\": <run-low>, \"dc\": <low-dc>},\n", "# 1, {\"gain\": 0, \"run_number\": <run-high>, \"dc\": <high-dc>},\n", "# 1, {\"gain\": 1, \"run_number\": <run-med>, \"dc\": <med-dc>},\n", "# 1, {\"gain\": 2, \"run_number\": <run-low>, \"dc\": <low-dc>},\n", "# ]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "offset_g = OrderedDict()\n", "noise_g = OrderedDict()\n", "badpix_g = OrderedDict()\n", "if not fixed_gain_mode:\n", " gain_g = OrderedDict()\n", " gainstd_g = OrderedDict()\n", "\n", "\n", "for module_index, gain_index, offset, noise, gains, gains_std, bp in results:\n", " qm = module_index_to_qm(module_index)\n", " if qm not in offset_g:\n", " offset_g[qm] = np.zeros((offset.shape[0], offset.shape[1], offset.shape[2], 3))\n", " noise_g[qm] = np.zeros_like(offset_g[qm])\n", " badpix_g[qm] = np.zeros_like(offset_g[qm], np.uint32)\n", " if not fixed_gain_mode:\n", " gain_g[qm] = np.zeros_like(offset_g[qm])\n", " gainstd_g[qm] = np.zeros_like(offset_g[qm])\n", "\n", " offset_g[qm][..., gain_index] = offset\n", " noise_g[qm][..., gain_index] = noise\n", " badpix_g[qm][..., gain_index] = bp\n", " if not fixed_gain_mode:\n", " gain_g[qm][..., gain_index] = gains\n", " gainstd_g[qm][..., gain_index] = gains_std" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Add bad pixels due to bad gain separation\n", "if not fixed_gain_mode:\n", " for qm in gain_g.keys():\n", " for g in range(2):\n", " # Bad pixels during bad gain separation.\n", " # Fraction of pixels in the module with separation lower than \"thresholds_gain_sigma\".\n", " bad_sep = (gain_g[qm][..., g+1] - gain_g[qm][..., g]) / \\\n", " np.sqrt(gainstd_g[qm][..., g+1]**2 + gainstd_g[qm][..., g]**2)\n", " badpix_g[qm][...,g+1][bad_sep<thresholds_gain_sigma] |= \\\n", " BadPixels.GAIN_THRESHOLDING_ERROR" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The thresholds for gain switching are then defined as the mean value between in individual gain bit levels. Note that these thresholds need to be refined with charge induced thresholds, as the two are not the same." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if not fixed_gain_mode:\n", " thresholds_g = {}\n", " for qm in gain_g.keys():\n", " thresholds_g[qm] = np.zeros((gain_g[qm].shape[0], gain_g[qm].shape[1], gain_g[qm].shape[2], 5))\n", " thresholds_g[qm][...,0] = (gain_g[qm][...,1]+gain_g[qm][...,0])/2\n", " thresholds_g[qm][...,1] = (gain_g[qm][...,2]+gain_g[qm][...,1])/2\n", " for i in range(3):\n", " thresholds_g[qm][...,2+i] = gain_g[qm][...,i]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res = OrderedDict()\n", "for i in modules:\n", " qm = module_index_to_qm(i)\n", " res[qm] = {\n", " 'Offset': offset_g[qm],\n", " 'Noise': noise_g[qm],\n", " 'BadPixelsDark': badpix_g[qm]\n", " }\n", " if not fixed_gain_mode:\n", " res[qm]['ThresholdsDark'] = thresholds_g[qm]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# set the operating condition\n", "# note: iCalibrationDB only adds gain_mode if it is truthy, so we don't need to handle None\n", "condition = iCalibrationDB.Conditions.Dark.AGIPD(\n", " memory_cells=mem_cells,\n", " bias_voltage=bias_voltage,\n", " acquisition_rate=acq_rate,\n", " gain_setting=gain_setting,\n", " gain_mode=fixed_gain_mode,\n", " integration_time=integration_time\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create mapping from module(s) (qm) to karabo_da(s) and PDU(s)\n", "qm_dict = OrderedDict()\n", "all_pdus = get_pdu_from_db(\n", " karabo_id,\n", " karabo_da,\n", " constant=iCalibrationDB.CalibrationConstant(),\n", " condition=condition,\n", " cal_db_interface=cal_db_interface,\n", " snapshot_at=creation_time.isoformat() if creation_time else None,\n", " timeout=cal_db_timeout\n", ")\n", "for module_index, module_da, module_pdu in zip(modules, karabo_da, all_pdus):\n", " qm = module_index_to_qm(module_index)\n", " qm_dict[qm] = {\n", " \"karabo_da\": module_da,\n", " \"db_module\": module_pdu\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sending calibration constants to the database." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "md = None\n", "\n", "for qm in res:\n", " db_module = qm_dict[qm][\"db_module\"]\n", " for const in res[qm]:\n", " dconst = getattr(iCalibrationDB.Constants.AGIPD, const)()\n", " dconst.data = res[qm][const]\n", "\n", " if db_output:\n", " md = send_to_db(db_module, karabo_id, dconst, condition, file_loc,\n", " report, cal_db_interface, 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, dconst.data,\n", " file_loc, report, creation_time, out_folder)\n", " print(f\"Calibration constant {const} for {qm} is stored locally in {file_loc}.\\n\")\n", "\n", "print(\"Constants parameter conditions are:\\n\")\n", "print(f\"• memory_cells: {mem_cells}\\n• bias_voltage: {bias_voltage}\\n\"\n", " f\"• acquisition_rate: {acq_rate}\\n• gain_setting: {gain_setting}\\n\"\n", " f\"• gain_mode: {fixed_gain_mode}\\n• integration_time: {integration_time}\\n\"\n", " f\"• creation_time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\\n\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Retrieving previous calibration constants for comparison." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Start retrieving existing constants for comparison\n", "qm_x_const = [(qm, const) for const in res[qm] for qm in res]\n", "\n", "\n", "def retrieve_old_constant(qm, const):\n", " dconst = getattr(iCalibrationDB.Constants.AGIPD, const)()\n", "\n", " data, mdata = get_from_db(\n", " karabo_id=karabo_id,\n", " karabo_da=qm_dict[qm][\"karabo_da\"],\n", " constant=dconst,\n", " condition=condition,\n", " empty_constant=None,\n", " cal_db_interface=cal_db_interface,\n", " creation_time=creation_time-timedelta(seconds=1) if creation_time else None,\n", " strategy=\"pdu_prior_in_time\",\n", " verbosity=1,\n", " timeout=cal_db_timeout\n", " )\n", "\n", " if mdata is None or data is None:\n", " timestamp = \"Not found\"\n", " filepath = None\n", " h5path = None\n", " else:\n", " timestamp = mdata.calibration_constant_version.begin_at.isoformat()\n", " filepath = os.path.join(\n", " mdata.calibration_constant_version.hdf5path,\n", " mdata.calibration_constant_version.filename\n", " )\n", " h5path = mdata.calibration_constant_version.h5path\n", "\n", " return data, timestamp, filepath, h5path\n", "\n", "\n", "old_retrieval_pool = multiprocessing.Pool()\n", "old_retrieval_res = old_retrieval_pool.starmap_async(\n", " retrieve_old_constant, qm_x_const\n", ")\n", "old_retrieval_pool.close()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mnames=[]\n", "for i in modules:\n", " qm = module_index_to_qm(i)\n", " mnames.append(qm)\n", " display(Markdown(f'## Position of the module {qm} and its ASICs'))\n", "show_processed_modules(dinstance, constants=None, mnames=mnames, mode=\"position\")" ] }, { "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 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": "markdown", "metadata": {}, "source": [ "### High Gain ###" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cell = 3\n", "gain = 0\n", "show_overview(res, cell, gain, infix=\"{}-{}-{}\".format(*run_numbers))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Medium Gain ###" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cell = 3\n", "gain = 1\n", "show_overview(res, cell, gain, infix=\"{}-{}-{}\".format(*run_numbers))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Low Gain ###" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cell = 3\n", "gain = 2\n", "show_overview(res, cell, gain, infix=\"{}-{}-{}\".format(*run_numbers))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if high_res_badpix_3d:\n", " cols = {\n", " BadPixels.NOISE_OUT_OF_THRESHOLD: (BadPixels.NOISE_OUT_OF_THRESHOLD.name, '#FF000080'),\n", " BadPixels.OFFSET_NOISE_EVAL_ERROR: (BadPixels.OFFSET_NOISE_EVAL_ERROR.name, '#0000FF80'),\n", " BadPixels.OFFSET_OUT_OF_THRESHOLD: (BadPixels.OFFSET_OUT_OF_THRESHOLD.name, '#00FF0080'),\n", " BadPixels.GAIN_THRESHOLDING_ERROR: (BadPixels.GAIN_THRESHOLDING_ERROR.name, '#FF40FF40'),\n", " BadPixels.OFFSET_OUT_OF_THRESHOLD | BadPixels.NOISE_OUT_OF_THRESHOLD: ('OFFSET_OUT_OF_THRESHOLD + NOISE_OUT_OF_THRESHOLD', '#DD00DD80'),\n", " BadPixels.OFFSET_OUT_OF_THRESHOLD | BadPixels.NOISE_OUT_OF_THRESHOLD |\n", " BadPixels.GAIN_THRESHOLDING_ERROR: ('MIXED', '#BFDF009F')\n", " }\n", "\n", " display(Markdown(\"\"\"\n", "\n", " ## Global Bad Pixel Behaviour ##\n", "\n", " The following plots show the results of bad pixel evaluation for all evaluated memory cells.\n", " Cells are stacked in the Z-dimension, while pixels values in x/y are rebinned 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", "\n", " gnames = ['High Gain', 'Medium Gain', 'Low Gain']\n", " for gain in range(3):\n", " display(Markdown(f'### {gnames[gain]} ###'))\n", " for mod, data in badpix_g.items():\n", " plot_badpix_3d(data[...,gain], cols, title=mod, rebin_fac=1)\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## 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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "create_constant_overview(offset_g, \"Offset (ADU)\", mem_cells, 4000, 8000,\n", " badpixels=[badpix_g, np.nan])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "create_constant_overview(noise_g, \"Noise (ADU)\", mem_cells, 0, 100,\n", " badpixels=[badpix_g, np.nan])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if not fixed_gain_mode:\n", " # Plot only three gain threshold maps.\n", " bp_thresh = OrderedDict()\n", " for mod, con in badpix_g.items():\n", " bp_thresh[mod] = np.zeros((con.shape[0], con.shape[1], con.shape[2], 5), dtype=con.dtype)\n", " bp_thresh[mod][...,:2] = con[...,:2]\n", " bp_thresh[mod][...,2:] = con\n", "\n", " create_constant_overview(thresholds_g, \"Threshold (ADU)\", mem_cells, 4000, 10000, 5,\n", " badpixels=[bp_thresh, np.nan],\n", " gmap=['HG-MG Threshold', 'MG-LG Threshold', 'High gain', 'Medium gain', 'low gain'],\n", " marker=['d','d','','','']\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bad_pixel_aggregate_g = OrderedDict()\n", "for m, d in badpix_g.items():\n", " bad_pixel_aggregate_g[m] = d.astype(np.bool).astype(np.float)\n", "create_constant_overview(bad_pixel_aggregate_g, \"Bad pixel fraction\", mem_cells, 0, 0.10, 3)" ] }, { "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": [ "# now we need the old constants\n", "old_const = {}\n", "old_mdata = {}\n", "old_retrieval_res.wait()\n", "\n", "for (qm, const), (data, timestamp, filepath, h5path) in zip(qm_x_const, old_retrieval_res.get()):\n", " old_const.setdefault(qm, {})[const] = data\n", " old_mdata.setdefault(qm, {})[const] = {\n", " \"timestamp\": timestamp,\n", " \"filepath\": filepath,\n", " \"h5path\": h5path\n", " }" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "display(Markdown(\"The following pre-existing constants are used for comparison:\"))\n", "for qm, consts in old_mdata.items():\n", " display(Markdown(f\"- {qm}\"))\n", " for const in consts:\n", " display(Markdown(f\" - {const} at {consts[const]['timestamp']}\"))\n", " # saving locations of old constants for summary notebook\n", " with open(f\"{out_folder}/module_metadata_{qm}.yml\", \"w\") as fd:\n", " yaml.safe_dump(\n", " {\n", " \"module\": qm,\n", " \"pdu\": qm_dict[qm][\"db_module\"],\n", " \"old-constants\": old_mdata[qm]\n", " },\n", " fd,\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "table = []\n", "gain_names = ['High', 'Medium', 'Low']\n", "bits = [BadPixels.NOISE_OUT_OF_THRESHOLD, BadPixels.OFFSET_OUT_OF_THRESHOLD, BadPixels.OFFSET_NOISE_EVAL_ERROR, BadPixels.GAIN_THRESHOLDING_ERROR]\n", "for qm in badpix_g.keys():\n", " for gain in range(3):\n", " l_data = []\n", " l_data_old = []\n", "\n", " data = np.copy(badpix_g[qm][:,:,:,gain])\n", " datau32 = data.astype(np.uint32)\n", " l_data.append(len(datau32[datau32>0].flatten()))\n", " for bit in bits:\n", " l_data.append(np.count_nonzero(badpix_g[qm][:,:,:,gain] & bit))\n", "\n", " if old_const[qm]['BadPixelsDark'] is not None:\n", " dataold = np.copy(old_const[qm]['BadPixelsDark'][:, :, :, gain])\n", " datau32old = dataold.astype(np.uint32)\n", " l_data_old.append(len(datau32old[datau32old>0].flatten()))\n", " for bit in bits:\n", " l_data_old.append(np.count_nonzero(old_const[qm]['BadPixelsDark'][:, :, :, gain] & bit))\n", "\n", " l_data_name = ['All bad pixels', 'NOISE_OUT_OF_THRESHOLD',\n", " 'OFFSET_OUT_OF_THRESHOLD', 'OFFSET_NOISE_EVAL_ERROR', 'GAIN_THRESHOLDING_ERROR']\n", "\n", " l_threshold = ['', f'{thresholds_noise_sigma}' f'{thresholds_noise_hard[gain]}',\n", " f'{thresholds_offset_sigma}' f'{thresholds_offset_hard[gain]}',\n", " '', f'{thresholds_gain_sigma}']\n", "\n", " for i in range(len(l_data)):\n", " line = [f'{l_data_name[i]}, {gain_names[gain]} 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", "### 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", " \"New constant\", \"Old constant \"]\n", "\n", "if fixed_gain_mode:\n", " constants = ['Offset', 'Noise']\n", "else:\n", " constants = ['Offset', 'Noise', 'ThresholdsDark']\n", "\n", "constants_x_qms = list(itertools.product(constants, res.keys()))\n", "\n", "\n", "def compute_table(const, qm):\n", " if const == 'ThresholdsDark':\n", " table = [['','HG-MG threshold', 'HG-MG threshold', 'MG-LG threshold', 'MG-LG threshold']]\n", " else:\n", " table = [['','High gain', 'High gain', 'Medium gain', 'Medium gain', 'Low gain', 'Low gain']]\n", "\n", " compare_with_old_constant = old_const[qm][const] is not None and \\\n", " old_const[qm]['BadPixelsDark'] is not None\n", "\n", " data = np.copy(res[qm][const])\n", "\n", " if const == 'ThresholdsDark':\n", " data[...,0][res[qm]['BadPixelsDark'][...,0]>0] = np.nan\n", " data[...,1][res[qm]['BadPixelsDark'][...,1]>0] = np.nan\n", " else:\n", " data[res[qm]['BadPixelsDark']>0] = np.nan\n", "\n", " if compare_with_old_constant:\n", " data_old = np.copy(old_const[qm][const])\n", " if const == 'ThresholdsDark':\n", " data_old[...,0][old_const[qm]['BadPixelsDark'][...,0]>0] = np.nan\n", " data_old[...,1][old_const[qm]['BadPixelsDark'][...,1]>0] = np.nan\n", " else:\n", " data_old[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", " def compute_row(i):\n", " line = [n_list[i]]\n", " for gain in range(3):\n", " # Compare only 3 threshold gain-maps\n", " if gain == 2 and const == 'ThresholdsDark':\n", " continue\n", " stat_measure = f_list[i](data[...,gain])\n", " line.append(f\"{stat_measure:6.1f}\")\n", " if compare_with_old_constant:\n", " old_stat_measure = f_list[i](data_old[...,gain])\n", " line.append(f\"{old_stat_measure:6.1f}\")\n", " else:\n", " line.append(\"-\")\n", " return line\n", "\n", "\n", " with multiprocessing.pool.ThreadPool(processes=multiprocessing.cpu_count() // len(constants_x_qms)) as pool:\n", " rows = pool.map(compute_row, range(len(f_list)))\n", "\n", " table.extend(rows)\n", "\n", " return table\n", "\n", "\n", "with multiprocessing.Pool(processes=len(constants_x_qms)) as pool:\n", " tables = pool.starmap(compute_table, constants_x_qms)\n", "\n", "for (const, qm), table in zip(constants_x_qms, tables):\n", " display(Markdown(f\"### {qm}: {const} [ADU], good pixels only\"))\n", " display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=header)))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "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.12" } }, "nbformat": 4, "nbformat_minor": 4 }