{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AGIPD Characterize Dark Images #\n",
    "\n",
    "Author: S. Hauf, Version: 0.1\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 = \"/gpfs/exfel/data/scratch/hammerd/agipd-fixed-gain\" # path to output to, required\n",
    "sequences = [0] # sequence files to evaluate.\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_id = \"{}CH0\" # inset for receiver devices\n",
    "path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data\n",
    "h5path = '/INSTRUMENT/{}/DET/{}:xtdf/image' # path in the HDF5 file to images\n",
    "h5path_idx = '/INDEX/{}/DET/{}:xtdf/image' # path in the HDF5 file to images\n",
    "h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP' # path to control information\n",
    "karabo_id_control = \"HED_EXP_AGIPD500K2G\" # karabo-id for control device '\n",
    "karabo_da_control = \"AGIPD500K2G00\" # karabo DA for control infromation\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 # detector bias voltage\n",
    "gain_setting = 0.1 # the gain setting, use 0.1 to try to auto-determine\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",
    "rawversion = 2 # RAW file format version\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",
    "\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "from datetime import datetime\n",
    "\n",
    "import dateutil.parser\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "import os\n",
    "from collections import OrderedDict\n",
    "from typing import List, Tuple\n",
    "\n",
    "import h5py\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import tabulate\n",
    "from cal_tools.enums import BadPixels\n",
    "\n",
    "matplotlib.use('agg')\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import Latex, Markdown, display\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import multiprocessing\n",
    "\n",
    "from cal_tools.agipdlib import (get_acq_rate, get_bias_voltage, get_gain_mode,\n",
    "                                get_gain_setting, get_num_cells)\n",
    "from cal_tools.enums import AgipdGainMode\n",
    "from cal_tools.plotting import (create_constant_overview, plot_badpix_3d,\n",
    "                                show_overview, show_processed_modules)\n",
    "from cal_tools.tools import (get_dir_creation_date, get_from_db,\n",
    "                             get_pdu_from_db, get_random_db_interface,\n",
    "                             get_report, map_gain_stages, module_index_to_qm,\n",
    "                             run_prop_seq_from_path, save_const_to_h5,\n",
    "                             send_to_db)\n",
    "from iCalibrationDB import Conditions, Constants, Detectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert control device if format string (does nothing otherwise)\n",
    "h5path_ctrl = h5path_ctrl.format(karabo_id_control)\n",
    "\n",
    "max_cells = mem_cells\n",
    "\n",
    "offset_runs = OrderedDict()\n",
    "offset_runs[\"high\"] = run_high\n",
    "offset_runs[\"med\"] = run_med\n",
    "offset_runs[\"low\"] = run_low\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",
    "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",
    "control_names = [f'{in_folder}/r{r:04d}/RAW-R{r:04d}-{karabo_da_control}-S00000.h5'\n",
    "                 for r in (run_high, run_med, run_low)]\n",
    "\n",
    "if operation_mode not in (\"ADAPTIVE_GAIN\", \"FIXED_GAIN\"):\n",
    "    print(f\"WARNING: unknown operation_mode \\\"{operation_mode}\\\" parameter set\")\n",
    "run_gain_modes = [get_gain_mode(fn, h5path_ctrl) for fn in control_names]\n",
    "if all(gm == AgipdGainMode.ADAPTIVE_GAIN for gm in run_gain_modes):\n",
    "    fixed_gain_mode = False\n",
    "    if operation_mode == \"FIXED_GAIN\":\n",
    "        print(\"WARNING: operation_mode parameter is FIXED_GAIN, slow data indicates adaptive gain\")\n",
    "elif run_gain_modes == [AgipdGainMode.FIXED_HIGH_GAIN, AgipdGainMode.FIXED_MEDIUM_GAIN, AgipdGainMode.FIXED_LOW_GAIN]:\n",
    "    if operation_mode == \"ADAPTIVE_GAIN\":\n",
    "        print(\"WARNING: operation_mode parameter ix ADAPTIVE_GAIN, slow data indicates fixed gain\")\n",
    "    fixed_gain_mode = True\n",
    "else:\n",
    "    print(f'Something is clearly wrong; slow data indicates gain modes {run_gain_modes}')\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": [
    "runs = [run_high, run_med, run_low]\n",
    "\n",
    "if gain_setting == 0.1:\n",
    "    if creation_time.replace(tzinfo=None) < dateutil.parser.parse('2020-01-31'):\n",
    "        print(\"Set gain-setting to None for runs taken before 2020-01-31\")\n",
    "        gain_setting = None\n",
    "    else:\n",
    "        try:\n",
    "            # extract gain setting and validate that all runs have the same setting\n",
    "            gsettings = []\n",
    "            for r in runs:\n",
    "                control_fname = '{}/r{:04d}/RAW-R{:04d}-{}-S00000.h5'.format(in_folder, r, r,\n",
    "                                                                             karabo_da_control)\n",
    "                gsettings.append(get_gain_setting(control_fname, h5path_ctrl))\n",
    "            if not all(g == gsettings[0] for g in gsettings):\n",
    "                raise ValueError(f\"Different gain settings for the 3 input runs {gsettings}\")\n",
    "            gain_setting =  gsettings[0]\n",
    "        except Exception as e:\n",
    "            print(f'Error while reading gain setting from: \\n{control_fname}')\n",
    "            print(f'Error: {e}')\n",
    "            if \"component not found\" in str(e):\n",
    "                print(\"Gain setting is not found in the control information\")\n",
    "            print(\"Data will not be processed\")\n",
    "            sequences = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if karabo_da[0] == '-1':\n",
    "    if modules[0] == -1:\n",
    "        modules = list(range(nmods))\n",
    "    karabo_da = [\"AGIPD{:02d}\".format(i) for i in modules]\n",
    "else:\n",
    "    modules = [int(x[-2:]) for x in karabo_da]\n",
    "h5path = h5path.format(karabo_id, receiver_id)\n",
    "h5path_idx = h5path_idx.format(karabo_id, receiver_id)\n",
    "\n",
    "if bias_voltage == 0:\n",
    "    # Read the bias voltage from files, if recorded.\n",
    "    # If not available, make use of the historical voltage the detector is running at\n",
    "    bias_voltage = get_bias_voltage(control_names[0], karabo_id_control)\n",
    "    bias_voltage = bias_voltage if bias_voltage is not None else 300\n",
    "\n",
    "print(\"Parameters are:\")\n",
    "print(f\"Proposal: {prop}\")\n",
    "print(f\"Memory cells: {mem_cells}/{max_cells}\")\n",
    "print(\"Runs: {}\".format([ v for v in offset_runs.values()]))\n",
    "print(f\"Sequences: {sequences}\")\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\"Operation mode is {'fixed' if fixed_gain_mode else 'adaptive'} gain mode\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following lines will create a queue of files which will the be executed module-parallel. Distiguishing between different gains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set everything up filewise\n",
    "os.makedirs(out_folder, exist_ok=True)\n",
    "gmf = map_gain_stages(in_folder, offset_runs, path_template, karabo_da, sequences)\n",
    "gain_mapped_files, total_sequences, total_file_size = gmf\n",
    "print(f\"Will process a total of {total_sequences} files.\")"
   ]
  },
  {
   "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": [
    "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(f\"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": [
    "def characterize_module(fast_data_filename: str, channel: int, gg: int) -> Tuple[np.array, np.array, np.array, np.array, int, np.array, int, float]:\n",
    "    if max_cells == 0:\n",
    "        num_cells = get_num_cells(fast_data_filename, karabo_id, channel)\n",
    "    else:\n",
    "        num_cells = max_cells\n",
    "\n",
    "    print(f\"Using {num_cells} memory cells\")\n",
    "\n",
    "    if acq_rate == 0.:\n",
    "        slow_paths = control_names[gg], karabo_id_control\n",
    "        fast_paths = fast_data_filename, karabo_id, channel\n",
    "        local_acq_rate = get_acq_rate(fast_paths, slow_paths)\n",
    "    else:\n",
    "        local_acq_rate = acq_rate\n",
    "\n",
    "    local_thresholds_offset_hard = thresholds_offset_hard[gg]\n",
    "    local_thresholds_noise_hard = thresholds_noise_hard[gg]\n",
    "\n",
    "    h5path_f = h5path.format(channel)\n",
    "    h5path_idx_f = h5path_idx.format(channel)\n",
    "\n",
    "    with h5py.File(fast_data_filename, \"r\") as infile:\n",
    "        if rawversion == 2:\n",
    "            count = np.squeeze(infile[f\"{h5path_idx_f}/count\"])\n",
    "            first = np.squeeze(infile[f\"{h5path_idx_f}/first\"])\n",
    "            last_index = int(first[count != 0][-1]+count[count != 0][-1])\n",
    "            first_index = int(first[count != 0][0])\n",
    "        else:\n",
    "            status = np.squeeze(infile[f\"{h5path_idx_f}/status\"])\n",
    "            if np.count_nonzero(status != 0) == 0:\n",
    "                return\n",
    "            last = np.squeeze(infile[f\"{h5path_idx_f}/last\"])\n",
    "            first = np.squeeze(infile[f\"{h5path_idx_f}/first\"])\n",
    "            last_index = int(last[status != 0][-1]) + 1\n",
    "            first_index = int(first[status != 0][0])\n",
    "        im = np.array(infile[f\"{h5path_f}/data\"][first_index:last_index,...])\n",
    "        cellIds = np.squeeze(infile[f\"{h5path_f}/cellId\"][first_index:last_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",
    "        cellIds = cellIds[::2]\n",
    "    else:\n",
    "        if not fixed_gain_mode:\n",
    "            ga = im[:, 1, ...]\n",
    "        im = im[:, 0, ...].astype(np.float32)\n",
    "\n",
    "    im = np.rollaxis(im, 2)\n",
    "    im = np.rollaxis(im, 2, 1)\n",
    "\n",
    "    if not fixed_gain_mode:\n",
    "        ga = np.rollaxis(ga, 2)\n",
    "        ga = np.rollaxis(ga, 2, 1)\n",
    "    \n",
    "    offset = np.zeros((im.shape[0], im.shape[1], num_cells))\n",
    "    noise = np.zeros((im.shape[0], im.shape[1], num_cells))\n",
    "\n",
    "    if fixed_gain_mode:\n",
    "        gains = None\n",
    "        gains_std = None\n",
    "    else:\n",
    "        gains = np.zeros((im.shape[0], im.shape[1], num_cells))\n",
    "        gains_std = np.zeros((im.shape[0], im.shape[1], num_cells))\n",
    "\n",
    "    for cc in np.unique(cellIds[cellIds < num_cells]):\n",
    "        cellidx = cellIds == cc\n",
    "        offset[...,cc] = np.median(im[..., cellidx], axis=2)\n",
    "        noise[...,cc] = np.std(im[..., cellidx], axis=2)\n",
    "        if not fixed_gain_mode:\n",
    "            gains[...,cc] = np.median(ga[..., cellidx], axis=2)\n",
    "            gains_std[...,cc] = np.std(ga[..., cellidx], axis=2)\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 < local_thresholds_offset_hard[0]) | (\n",
    "        offset > local_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",
    "    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 < local_thresholds_noise_hard[0]) | (noise > local_thresholds_noise_hard[1])] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value\n",
    "    bp[~np.isfinite(noise)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value\n",
    "\n",
    "    return offset, noise, gains, gains_std, gg, bp, num_cells, local_acq_rate"
   ]
  },
  {
   "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",
    "start = datetime.now()\n",
    "all_cells = []\n",
    "all_acq_rate = []\n",
    "\n",
    "inp = []\n",
    "for gg, (gain, mapped_files) in enumerate(gain_mapped_files.items()):\n",
    "    dones = []\n",
    "    for i in modules:\n",
    "        qm = module_index_to_qm(i)\n",
    "        if qm in mapped_files and not mapped_files[qm].empty():\n",
    "            fname_in = mapped_files[qm].get()\n",
    "            print(\"Process file: \", fname_in)\n",
    "            dones.append(mapped_files[qm].empty())\n",
    "        else:\n",
    "            continue\n",
    "        inp.append((fname_in, i, gg))\n",
    "\n",
    "with multiprocessing.Pool() as pool:\n",
    "    results = pool.starmap(characterize_module, inp)\n",
    "\n",
    "for offset, noise, gains, gains_std, gg, bp, thiscell, thisacq in results:\n",
    "    all_cells.append(thiscell)\n",
    "    all_acq_rate.append(thisacq)\n",
    "    for i in modules:\n",
    "        qm = module_index_to_qm(i)\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][...,gg] = offset\n",
    "        noise_g[qm][...,gg] = noise\n",
    "        badpix_g[qm][...,gg] = bp\n",
    "        if not fixed_gain_mode:\n",
    "            gain_g[qm][...,gg] = gains\n",
    "            gainstd_g[qm][..., gg] = gains_std\n",
    "\n",
    "\n",
    "duration = (datetime.now() - start).total_seconds()\n",
    "\n",
    "max_cells = np.max(all_cells)\n",
    "print(f\"Using {max_cells} memory cells\")\n",
    "acq_rate = np.max(all_acq_rate)\n",
    "print(f\"Using {acq_rate} MHz acquisition rate\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add a badpixel due to bad gain separation\n",
    "if not fixed_gain_mode:\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]) / np.sqrt(gainstd_g[qm][..., g+1]**2 + gainstd_g[qm][..., g]**2)\n",
    "        badpix_g[qm][...,g+1][(bad_sep)<thresholds_gain_sigma]|= BadPixels.GAIN_THRESHOLDING_ERROR.value"
   ]
  },
  {
   "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": [
    "# 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(out_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 = OrderedDict()\n",
    "for i, k_da in zip(modules, karabo_da):\n",
    "    qm = module_index_to_qm(i)\n",
    "    qm_dict[qm] = {\n",
    "        \"karabo_da\": k_da,\n",
    "        \"db_module\": \"\"\n",
    "    }"
   ]
  },
  {
   "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 = Conditions.Dark.AGIPD(memory_cells=max_cells,\n",
    "                                  bias_voltage=bias_voltage,\n",
    "                                  acquisition_rate=acq_rate,\n",
    "                                  gain_setting=gain_setting,\n",
    "                                  gain_mode=fixed_gain_mode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve existing constants for comparison\n",
    "if fixed_gain_mode:\n",
    "    clist = [\"Offset\", \"Noise\", \"BadPixelsDark\"]\n",
    "else:\n",
    "    clist = [\"Offset\", \"Noise\", \"ThresholdsDark\", \"BadPixelsDark\"]\n",
    "old_const = {}\n",
    "old_mdata = {}\n",
    "detinst = getattr(Detectors, dinstance)\n",
    "\n",
    "print('Retrieve pre-existing constants for comparison.')\n",
    "for qm in res:\n",
    "    qm_db = qm_dict[qm]\n",
    "    karabo_da = qm_db[\"karabo_da\"]\n",
    "    for const in res[qm]:\n",
    "        dconst = getattr(Constants.AGIPD, const)()\n",
    "        dconst.data = res[qm][const]\n",
    "\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",
    "        # TODO: Set db_module to \"\" by default in the first cell\n",
    "        if not qm_db[\"db_module\"]:\n",
    "            qm_db[\"db_module\"] = get_pdu_from_db(karabo_id, karabo_da, dconst,\n",
    "                                                 condition, cal_db_interface,\n",
    "                                                 snapshot_at=creation_time)[0]\n",
    "\n",
    "        data, mdata = get_from_db(karabo_id, karabo_da,\n",
    "                                  dconst, condition,\n",
    "                                  None, cal_db_interface,\n",
    "                                  creation_time=creation_time,\n",
    "                                  verbosity=2, timeout=cal_db_timeout)\n",
    "\n",
    "        old_const[const] = data\n",
    "        if mdata is not None and data is not None:\n",
    "            time = mdata.calibration_constant_version.begin_at\n",
    "            old_mdata[const] = time.isoformat()\n",
    "            os.makedirs('{}/old/'.format(out_folder), exist_ok=True)\n",
    "            save_const_to_h5(qm_db[\"db_module\"], karabo_id,\n",
    "                             dconst, condition, data,\n",
    "                             file_loc, report, creation_time,\n",
    "                             f'{out_folder}/old/')\n",
    "        else:\n",
    "            old_mdata[const] = \"Not found\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "md = None\n",
    "\n",
    "for qm in res:\n",
    "    karabo_da = qm_dict[qm][\"karabo_da\"]\n",
    "    db_module = qm_dict[qm][\"db_module\"]\n",
    "    for const in res[qm]:\n",
    "        dconst = getattr(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} is stored locally.\\n\")\n",
    "\n",
    "    print(\"Constants parameter conditions are:\\n\")\n",
    "    print(f\"• memory_cells: {max_cells}\\n• bias_voltage: {bias_voltage}\\n\"\n",
    "          f\"• acquisition_rate: {acq_rate}\\n• gain_setting: {gain_setting}\\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": [
    "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(*offset_runs.values()))"
   ]
  },
  {
   "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(*offset_runs.values()))"
   ]
  },
  {
   "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(*offset_runs.values()))"
   ]
  },
  {
   "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.GAIN_THRESHOLDING_ERROR.value: (BadPixels.GAIN_THRESHOLDING_ERROR.name, '#FF40FF40'),\n",
    "        BadPixels.OFFSET_OUT_OF_THRESHOLD.value | BadPixels.NOISE_OUT_OF_THRESHOLD.value: ('OFFSET_OUT_OF_THRESHOLD + NOISE_OUT_OF_THRESHOLD', '#DD00DD80'),\n",
    "        BadPixels.OFFSET_OUT_OF_THRESHOLD.value | BadPixels.NOISE_OUT_OF_THRESHOLD.value |\n",
    "        BadPixels.GAIN_THRESHOLDING_ERROR.value: ('MIXED', '#BFDF009F')}\n",
    "\n",
    "if high_res_badpix_3d:\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": [
    "## 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)\", max_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)\", max_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",
    "\n",
    "    create_constant_overview(thresholds_g, \"Threshold (ADU)\", max_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\", max_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": [
    "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",
    "\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.value))\n",
    "\n",
    "        if old_const['BadPixelsDark'] is not None:\n",
    "            dataold = np.copy(old_const['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['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', '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['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",
    "          \"New constant\", \"Old constant \"]\n",
    "\n",
    "if fixed_gain_mode:\n",
    "    constants = ['Offset', 'Noise']\n",
    "else:\n",
    "    constants = ['Offset', 'Noise', 'ThresholdsDark']\n",
    "\n",
    "for const in constants:\n",
    "\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",
    "    for qm in res.keys():\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 old_const[const] is not None and old_const['BadPixelsDark'] is not None:\n",
    "            dataold = np.copy(old_const[const])\n",
    "            if const == 'ThresholdsDark':\n",
    "                dataold[...,0][old_const['BadPixelsDark'][...,0]>0] = np.nan\n",
    "                dataold[...,1][old_const['BadPixelsDark'][...,1]>0] = np.nan\n",
    "            else:\n",
    "                dataold[old_const['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",
    "                # Compare only 3 threshold gain-maps\n",
    "                if gain == 2 and const == 'ThresholdsDark':\n",
    "                    continue\n",
    "                line.append('{:6.1f}'.format(f(data[...,gain])))\n",
    "                if old_const[const] is not None and old_const['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": "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}