diff --git a/notebooks/AGIPD/CS_Characterization_unequalClockStep_NBC.ipynb b/notebooks/AGIPD/CS_Characterization_unequalClockStep_NBC.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..53443e4c53182f49cc0237355ec85a942808ae8b --- /dev/null +++ b/notebooks/AGIPD/CS_Characterization_unequalClockStep_NBC.ipynb @@ -0,0 +1,1516 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Characterize AGIPD Current Source Data\n", + "\n", + "Author: Detector Group, Version 1.0\n", + "\n", + "The following notebook characterises the AGIPD data taken with the current source (CS) and is used to produce constants to allow for relative gain determination." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "in_folder='/gpfs/exfel/exp/SPB/202231/p900295/scratch/CSmergedFiles/07102022/' # path to input data, required\n", + "out_folder = \"/gpfs/exfel/exp/SPB/202231/p900295/scratch/CSconsts/14122022\" # path to output to, required\n", + "raw_folder = '/gpfs/exfel/exp/SPB/202231/p900295/raw/' # path to raw folder, required\n", + "first_run = 19\n", + "dark_run = 13\n", + "\n", + "modules = [4] # modules to work on, required, range allowed\n", + "karabo_da = [\"all\"]\n", + "karabo_da_control = \"AGIPD1MCTRL00\" # karabo DA for control infromation\n", + "karabo_id_control = \"SPB_IRU_AGIPD1M1\" # karabo-id for the control device e.g. \"MID_EXP_AGIPD1M1\", or \"SPB_IRU_AGIPD1M1\"\n", + "karabo_id = \"SPB_DET_AGIPD1M-1\"\n", + "ctrl_source_template = '{}/MDL/FPGA_COMP' # path to control information\n", + "instrument_source_template = '{}/DET/{}:xtdf' # path in the HDF5 file to images\n", + "receiver_template = \"{}CH0\" # inset for receiver devices\n", + "\n", + "use_dir_creation_date = True\n", + "delta_time = 0 # offset to the creation time (e.g. useful in case we want to force the system to use diff. dark constants)\n", + "cal_db_interface = \"tcp://max-exfl016:8019\" # the database interface to use\n", + "local_output = True # output constants locally\n", + "db_output = False # output constants to database\n", + "\n", + "bias_voltage = -1 # detector bias voltage\n", + "mem_cells = -1 # number of memory cells used, use -1 to auto-derive\n", + "acq_rate = 4.5 # the detector acquisition rate, use -1. to try to auto-determine\n", + "gain_setting = 0 # gain setting can have value 0 or 1, Default=0.1 for no (None) gain-setting\n", + "integration_time = 12 # integration time, negative values for auto-detection.\n", + "\n", + "sigma_dev_cut = 3 # parameters outside the range median +- sigma_dev_cut*MAD are replaced with the median\n", + "chi2_lim = 7 # limit of chi2 of the fit\n", + "relgain_lim = [0.7, 1.3] # limit for the relative gain\n", + "\n", + "steps = [1, 10, 50]\n", + "increments = [300, 500, 100]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import warnings\n", + "from datetime import datetime, timedelta\n", + "from functools import partial\n", + "\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "import dateutil.parser\n", + "import h5py\n", + "import matplotlib\n", + "import numpy as np\n", + "\n", + "from scipy.stats import median_abs_deviation as mad\n", + "from collections import OrderedDict\n", + "import xarray\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import gridspec\n", + "from matplotlib.colors import LogNorm, PowerNorm\n", + "import matplotlib.patches as patches\n", + "from mpl_toolkits.axes_grid1 import AxesGrid\n", + "\n", + "import XFELDetAna.xfelpyanatools as xana\n", + "from extra_data import RunDirectory, components\n", + "import pasha as psh\n", + "import multiprocessing\n", + "import itertools\n", + "from sklearn import mixture, cluster, preprocessing\n", + "\n", + "from cal_tools.agipdlib import AgipdCtrl\n", + "from cal_tools.enums import BadPixels\n", + "from cal_tools.plotting import plot_badpix_3d, show_overview\n", + "from cal_tools.tools import (\n", + " gain_map_files,\n", + " get_constant_from_db_and_time,\n", + " get_dir_creation_date,\n", + " get_notebook_name,\n", + " get_pdu_from_db,\n", + " get_report,\n", + " module_index_to_qm,\n", + " parse_runs,\n", + " send_to_db,\n", + ")\n", + "from iCalibrationDB import Conditions, ConstantMetaData, Constants, Detectors, Versions\n", + "\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for module in modules:\n", + " base = 'agipd_CH{:02d}'.format(module)\n", + " for fname in os.listdir(in_folder): \n", + " if base in fname:\n", + " print(f'Processing file: {fname}')\n", + " merged_h5file = fname\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cells = mem_cells\n", + "path_temp = raw_folder+\"/r{:04d}/\"\n", + "image_name_temp = 'RAW-R{:04d}-AGIPD{:02d}-S{:05d}.h5'\n", + "print(\"Parameters are:\")\n", + "if mem_cells < 0:\n", + " print(\"Memory cells: auto-detection on\")\n", + "else:\n", + " print(\"Memory cells set by user: {}\".format(mem_cells))\n", + "print(\"Run: {}\".format(first_run))\n", + "print(\"Modules: {}\".format(modules))\n", + "\n", + "instrument = karabo_id.split(\"_\")[0]\n", + "\n", + "\n", + "if instrument == \"HED\":\n", + " nmods = 8\n", + "else:\n", + " nmods = 16\n", + "\n", + "print(f\"Detector in use is {karabo_id}\")\n", + "\n", + "\n", + "if karabo_da == [\"all\"]:\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]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "channel = 0\n", + "ctrl_src = ctrl_source_template.format(karabo_id_control)\n", + "instrument_src = instrument_source_template.format(karabo_id, receiver_template).format(channel)\n", + "\n", + "agipd_cond = AgipdCtrl(\n", + " run_dc=RunDirectory(f'{raw_folder}/r{first_run:04d}'),\n", + " image_src=instrument_src,\n", + " ctrl_src=ctrl_src,\n", + " raise_error=False, # to be able to process very old data without gain_setting value\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define creation time\n", + "creation_time = None\n", + "if use_dir_creation_date:\n", + " creation_time = get_dir_creation_date(raw_folder, first_run)\n", + " creation_time = creation_time + timedelta(hours=delta_time)\n", + "\n", + "# Read AGIPD parameter conditions.\n", + "if acq_rate < 0.:\n", + " acq_rate = agipd_cond.get_acq_rate()\n", + "if mem_cells < 0:\n", + " mem_cells = agipd_cond.get_num_cells()\n", + "# TODO: look for alternative for passing creation_time\n", + "if gain_setting < 0:\n", + " gain_setting = agipd_cond.get_gain_setting(creation_time)\n", + "if bias_voltage < 0:\n", + " bias_voltage = agipd_cond.get_bias_voltage(karabo_id_control)\n", + "if integration_time < 0:\n", + " integration_time = agipd_cond.get_integration_time()\n", + "\n", + "print(f\"Acquisition rate: {acq_rate} MHz\")\n", + "print(f\"Memory cells: {mem_cells}\")\n", + "print(f\"Gain setting: {gain_setting}\")\n", + "print(f\"Integration time: {integration_time}\")\n", + "print(f\"Using {creation_time} as creation time of constant.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cs_data = []\n", + "with h5py.File(in_folder+merged_h5file, 'r') as f:\n", + " cs_data.append(np.array(f[f'/{modules[0]}/Analog/data']))\n", + "\n", + "trains = cs_data[0].shape[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a = np.arange(0,increments[0],steps[0])\n", + "b = np.arange(a[-1]+1,a[-1]+1+increments[1]*steps[1], steps[1])\n", + "c = np.arange(b[-1]+1,b[-1]+1+increments[2]*steps[2], steps[2])\n", + "x_eq = np.concatenate((a,b,c))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from iminuit import Minuit\n", + "from iminuit.util import describe, make_func_code\n", + "\n", + "def rolling_window(a, window):\n", + " shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)\n", + " strides = a.strides + (a.strides[-1],)\n", + " return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)\n", + "\n", + "def fit_data(fun, x, y, par_ests, minos_errs=False):\n", + " par_ests[\"throw_nan\"] = False\n", + " par_ests[\"pedantic\"] = False\n", + "# par_ests[\"print_level\"] = 0\n", + "\n", + " f_sig = describe(fun)[1:]\n", + "\n", + " class _Chi2Functor:\n", + " def __init__(self, f, x, y, err):\n", + " self.f = f\n", + " self.x = x[y != 0]\n", + " self.y = y[y != 0]\n", + " self.err = err[y != 0]\n", + " f_sig = describe(f)\n", + " # this is how you fake function\n", + " # signature dynamically\n", + " self.func_code = make_func_code(\n", + " f_sig[1:]) # docking off independent variable\n", + " self.func_defaults = None # this keeps numpy.vectorize happy\n", + "\n", + " def __call__(self, *arg):\n", + " # notice that it accept variable length\n", + " # positional arguments\n", + " # chi2 = sum((y-self.f(x,*arg))**2 for x,y in zip(self.x,self.y))\n", + " return np.sum(((self.f(self.x, *arg) - self.y) ** 2) / self.err)\n", + " \n", + " yerr = np.copy(y) #errors --> this has to be changed to use noise values instead\n", + " wrapped = _Chi2Functor(fun, x, y, yerr)\n", + " m = Minuit(wrapped, **par_ests)\n", + " m.migrad()\n", + " \n", + "\n", + " if m.get_fmin().is_valid:\n", + " if minos_errs:\n", + " m.minos()\n", + " print('{}\\n'.format(m.get_param_states()))\n", + " redchi2 = m.fval / (len(x) - len(m.args)) \n", + " else:\n", + " redchi2 = -1\n", + "# m.print_matrix() \n", + " fit_res = m.fitarg\n", + " fit_res['red_chi2'] = redchi2\n", + " \n", + " return fit_res\n", + "\n", + "def lin_fun(x, m, b):\n", + " return m*x+b" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Mark gain stage\n", + "Function \"label_gain_stage\" labels gain stage of provided data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def moving_average(vals, w):\n", + " return np.convolve(vals, np.ones(w), 'same') / w\n", + "\n", + "def label_gain_stage(x, y):\n", + " labels=[]\n", + " bounds = []\n", + " gain_marker = [0, 1, 2]\n", + " \n", + " if np.max(y) < 8e3:\n", + " pass\n", + " else:\n", + " diff = np.diff(y, prepend=y[0])\n", + "\n", + " mins = []\n", + " boundaries = []\n", + " mins.append(np.min(diff[x<300]))\n", + " mins.append(np.min(diff[x>300]))\n", + "\n", + " mins = np.array(mins)\n", + " if np.any(mins) > -200.:\n", + " mins = mins[mins < -200.]\n", + "\n", + " for m in mins:\n", + " boundaries.append(np.where(diff == m)[0][0])\n", + "\n", + " if (boundaries[-1] > 890) | (np.max(y) < 9e3): #| (np.max(y) > 14e3):\n", + " pass\n", + " elif len(mins) != 2:\n", + " pass\n", + " else:\n", + " hg_b = boundaries[0]\n", + " mg_b = boundaries[-1]\n", + "\n", + " labels = np.zeros(y.shape[0])\n", + " labels[:hg_b+1] = 0\n", + " labels[hg_b+1:mg_b+1] = 1\n", + " labels[mg_b+1:] = 2\n", + "\n", + " # Excluding boundary values shifted by 5 from each side\n", + " boundaries = np.array([np.arange(hg_b-5, hg_b+20), \n", + " np.arange(mg_b-5, mg_b+20), \n", + " ])\n", + " labels[boundaries] = -1\n", + " maxval_idx = np.where(y[labels==2] == np.max(y[labels==2]))[0][-1]\n", + " #######################################################\n", + " if (len(np.unique(labels)) < 4) | (np.min(y[labels==2]) > 13.5e3) | (maxval_idx < 50):\n", + " pass\n", + " else:\n", + " grad2 = np.gradient(y[labels==2])\n", + " mvavg = moving_average(grad2, 10)\n", + " grad = np.gradient(y[labels==1])\n", + "\n", + "\n", + " if len(mvavg) < 100:\n", + " msk2 = np.ones_like(mvavg, dtype=bool)\n", + " \n", + " else:\n", + " msk2 = 50+np.where(mvavg[50:maxval_idx]==np.max(mvavg[50:maxval_idx]))[0][-1]\n", + "\n", + "\n", + " msk = np.where(grad==np.max(grad[y[labels==1] > \n", + " 0.95*np.max(y[labels==1])]))[0][-1]\n", + " bounds.append(hg_b-6)\n", + " bounds.append(np.where(y[labels==1] == y[labels==1][msk])[0][0])\n", + " bounds.append(np.where(y[labels==2] == y[labels==2][msk2])[0][-1])\n", + "\n", + " return labels, gain_marker, bounds" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "run = RunDirectory(raw_folder+\"/r{:04d}/\".format(dark_run))\n", + "instrument_source = karabo_id+'/DET/{}CH0:xtdf'.format(modules[0])\n", + "offset = run.get_array(instrument_source, 'image.data')[:trains*mem_cells, 0, :, :]\n", + "offset = offset.values.reshape((offset.shape[0] // mem_cells, mem_cells, 1, 512, 128))\n", + "offset = np.moveaxis(offset, 4, 3)\n", + "offset = np.squeeze(offset).astype(np.float32)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# from cal_tools.tools import get_constant_from_db_and_time\n", + "\n", + "# noises = {}\n", + "\n", + "# for k_da, mod in zip(karabo_da, modules): \n", + "# noise, when = get_constant_from_db_and_time(karabo_id, k_da,\n", + "# Constants.AGIPD.Noise(),\n", + "# Conditions.Dark.AGIPD(\n", + "# memory_cells=mem_cells,\n", + "# bias_voltage=bias_voltage,\n", + "# acquisition_rate=acq_rate,\n", + "# gain_setting=gain_setting,\n", + "# integration_time=integration_time),\n", + "# np.zeros((128, 512, mem_cells, 3)),\n", + "# cal_db_interface, creation_time=creation_time)\n", + "# noises[mod] = np.array(noise.data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# %%time\n", + "test_pixels = []\n", + "tpix_range1 = [(200,210), (60,73)]\n", + "for i in range(*tpix_range1[0]):\n", + " for j in range(*tpix_range1[1]):\n", + " test_pixels.append((j,i))\n", + "test_cells = [1, 35]\n", + "tcell = np.array(test_cells)\n", + "tcell = tcell[tcell < mem_cells]\n", + "if tcell.size == 0:\n", + " test_cells = [mem_cells-1]\n", + "else:\n", + " test_cells = tcell.tolist()\n", + "\n", + "bins = (np.linspace(-100, 6000, 100),\n", + " np.linspace(5000, 13000, 1000),\n", + " ) \n", + " \n", + "from mpl_toolkits.axes_grid1 import ImageGrid\n", + "\n", + "ana = cs_data[0] \n", + "\n", + "H = [0, 0, 0, 0, 0]\n", + "gains = [0, 1, 1, 2, 2]\n", + "ex, ey = None, None\n", + "for pix in test_pixels:\n", + " for cell in test_cells:\n", + " color = np.random.rand(3,1)\n", + " x = x_eq[:ana.shape[0]]\n", + " y = ana[:,cell, pix[0], pix[1]]\n", + "\n", + " if np.any(y):\n", + " labels, gain_marker, bounds = label_gain_stage(x, y)\n", + "\n", + " if np.any(labels):\n", + " for i, lbl in enumerate(gain_marker):\n", + "\n", + " ym = y[labels==lbl]\n", + " h, ex, ey = np.histogram2d(x[labels==lbl], ym, range=((0,6e3), (5e3, 13e3)), bins=bins)\n", + " H[i] += h\n", + "\n", + "\n", + "fig = plt.figure(figsize=(10,7))\n", + "ax = fig.add_subplot(111)\n", + "ax.grid(lw=1.5)\n", + "for i in range(3):\n", + " H[i][H[i]==0] = np.nan \n", + "ax.imshow(H[0].T, origin=\"lower\", extent=[ex[0], ex[-1], ey[0], ey[-1]],\n", + " aspect='auto', cmap='summer', alpha=0.5, vmin=0, vmax=1000)\n", + "ax.imshow(H[1].T, origin=\"lower\", extent=[ex[0], ex[-1], ey[0], ey[-1]],\n", + " aspect='auto', cmap='spring', alpha=0.5, vmin=0, vmax=100)\n", + "ax.imshow(H[2].T, origin=\"lower\", extent=[ex[0], ex[-1], ey[0], ey[-1]],\n", + " aspect='auto', cmap='winter', alpha=0.5, vmin=0, vmax=1000)\n", + "ax.set_ylabel(\"AGIPD response (ADU)\", fontsize=12)\n", + "ax.set_xlabel(\"CS scan point (#)\", fontsize=12)\n", + "plt.xticks(fontsize=12)\n", + "plt.yticks(fontsize=12)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def evaluate_fitting_ROI(data, tpix_range, test_cells, roi):\n", + " lines1=0\n", + " no_entry = 0\n", + " fit_errors = {}\n", + " fits = {}\n", + " fits_wo = {}\n", + " reldiff = []\n", + " \n", + " markers = ['o', '.', 'x', 'v']\n", + " colors = ['tab:blue', 'tab:orange', 'tab:green']\n", + " test_pixels = []\n", + " \n", + " fig1, ax1 = plt.subplots(figsize=(9, 5))\n", + " ax1.grid(zorder=0, lw=1.5)\n", + " ax1.set_ylabel(\"CS signal (ADU)\", fontsize=11)\n", + " ax1.set_xlabel('Integration time (clk)', fontsize=11)\n", + " ax1.title.set_text('Raw CS signal')\n", + " \n", + " fig2, ax2 = plt.subplots(figsize=(9, 5))\n", + " ax2.grid(zorder=0, lw=1.5)\n", + " ax2.set_ylabel(\"CS signal (ADU) fit cut\", fontsize=11)\n", + " ax2.set_xlabel('Integration time (clk)', fontsize=11)\n", + " ax2.title.set_text('Offset corrected CS signal with fit cuts')\n", + " \n", + " fig3 = plt.figure(figsize=(9, 4))\n", + " gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1]) \n", + " ax3 = plt.subplot(gs[0])\n", + " ax4 = plt.subplot(gs[1])\n", + " ax3.grid(zorder=0, lw=1.5)\n", + " ax3.set_ylabel('CS signal (ADU)', fontsize=11)\n", + " ax4.set_ylabel('Relative\\ndeviation', fontsize=11)\n", + " ax4.set_xlabel('Integration time (clk)', fontsize=11)\n", + " ax3.title.set_text('CS signal fits')\n", + " \n", + " fig5 = plt.figure(figsize=(9, 4))\n", + " gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])\n", + " ax5 = plt.subplot(gs[0])\n", + " ax6 = plt.subplot(gs[1])\n", + " ax5.grid(zorder=0, lw=1.5)\n", + " ax6.set_xlabel(\"Pixel counter\", fontsize=11)\n", + " ax6.set_ylabel(\"Relative\\ndifference\", fontsize=11)\n", + " ax5.set_ylabel('Fitted slope ratio', fontsize=11)\n", + " ax5.title.set_text('Ratios of fitted gains')\n", + " ax5.set_ylim(0,60)\n", + " ax5.minorticks_on()\n", + " \n", + " fig7, ax7 = plt.subplots(figsize=(9, 5))\n", + " ax7.grid(zorder=0, lw=1.5)\n", + " ax7.set_xlabel('Integration time (clk)', fontsize=11)\n", + " ax7.set_ylabel('CS signal (keV)', fontsize=11)\n", + " ax7.title.set_text('Corrected CS signal (keV)')\n", + " \n", + " for i in range(*tpix_range[0]):\n", + " for j in range(*tpix_range[1]):\n", + " test_pixels.append((j,i))\n", + " \n", + " tcell = np.array(test_cells)\n", + " tcell = tcell[tcell < mem_cells]\n", + " if tcell.size == 0:\n", + " test_cells = [mem_cells-1]\n", + " else:\n", + " test_cells = tcell.tolist()\n", + " \n", + "\n", + " ana = data\n", + "# noise = noises[mod]\n", + " counter = 0\n", + " label_added = False\n", + "\n", + " for cell in test_cells:\n", + " for pix in test_pixels:\n", + " color = np.random.rand(3, 1)\n", + " x = x_eq[:ana.shape[0]]\n", + " y = ana[:trains,cell, pix[0], pix[1]]\n", + " \n", + " y_corr = (y - offset[:, cell, pix[0], pix[1]]).astype(np.float32) \n", + "\n", + " labels, gain_marker, bounds = label_gain_stage(x, y)\n", + " \n", + " if np.any(bounds):\n", + " for g in gain_marker:\n", + "\n", + " ax1.plot(x[labels==g][::4], y[labels==g][::4], \n", + " ls='None', marker=markers[g], color=colors[g], alpha=0.3)\n", + " nonzero = np.where(y_corr[labels==g] > 0.)\n", + " x_sg = x[labels==g][nonzero][:bounds[g]]\n", + " y_sg = y_corr[labels==g][nonzero][:bounds[g]]\n", + " y_wo = y[labels==g][nonzero][:bounds[g]]\n", + "\n", + " ax2.plot(x_sg, y_sg, ls='None', marker=markers[g], \n", + " color=colors[g], alpha=0.3)\n", + "\n", + " if y_sg.shape[0] != 0: \n", + " parms = {'m': 1, 'b': np.min(y_sg)}\n", + "\n", + "# errors = np.ones(x_sg.shape) * noise[pix[0], pix[1], cell, g]\n", + " fitted = fit_data(lin_fun, x_sg, y_sg, parms)\n", + " fitted_wo = fit_data(lin_fun, x_sg, y_wo, parms)\n", + " fits[g] = fitted\n", + " fits_wo[g] = fitted_wo\n", + "\n", + " yf = lin_fun(x_sg, fitted['m'], fitted['b'])\n", + "\n", + " ax3.plot(x_sg, y_sg, color=colors[g], ls='None',alpha=0.3, marker='o')\n", + " ax3.plot(x_sg, yf, color='navy', lw=1, zorder=9, alpha=0.2)\n", + " ax4.plot(x_sg, ((y_sg-yf) / y_sg), lw=1, color='navy')\n", + " hm = fits[0]['m'] / fits[1]['m']\n", + " hl = fits[0]['m'] / fits[2]['m']\n", + " ml = fits[1]['m'] / fits[2]['m']\n", + "\n", + " hm_wo = fits_wo[0]['m'] / fits_wo[1]['m']\n", + " ml_wo = fits_wo[1]['m'] / fits_wo[2]['m']\n", + "\n", + " M_offset = fits[1]['b']*hm-fits[0]['b']\n", + " L_offset = fits[2]['b']*hl-fits[0]['b']\n", + " y_mgcorr = y_corr[labels==1] * hm - M_offset\n", + " y_lgcorr = y_corr[labels==2] * hl - L_offset\n", + "\n", + " err_hm = hm * np.sqrt((fits[0]['error_m']/fits[0]['m'])**2 + (fits[1]['error_m']/fits[1]['m'])**2)\n", + " err_ml = ml * np.sqrt((fits[1]['error_m']/fits[1]['m'])**2 + (fits[2]['error_m']/fits[2]['m'])**2)\n", + "\n", + " if not label_added:\n", + " ax5.errorbar(counter, hm, err_hm, label='HG/MG', ls=\"None\", marker='x', color='tab:blue')\n", + " ax5.errorbar(counter, ml, err_ml, label='MG/LG', ls=\"None\", marker='x', color='tab:orange')\n", + " ax5.scatter(counter, hm_wo, label='w/o offset corr.', ls=\"None\", marker='x', color='tab:green')\n", + " ax5.scatter(counter, ml_wo, label='w/o offset corr.', ls=\"None\", marker='x', color='tab:red')\n", + " label_added = True\n", + " else:\n", + " ax5.errorbar(counter, hm, err_hm, ls=\"None\", marker='x', color='tab:blue')\n", + " ax5.errorbar(counter, ml, err_ml, ls=\"None\", marker='x', color='tab:orange')\n", + " ax5.scatter(counter, hm_wo, ls=\"None\", marker='x', color='tab:green')\n", + " ax5.scatter(counter, ml_wo, ls=\"None\", marker='x', color='tab:red')\n", + " reldiff.append((hm_wo - hm) / hm)\n", + " ax6.scatter(counter, ((hm_wo - hm) / hm), marker='.', color='tab:blue')\n", + "\n", + " counter += 1\n", + "\n", + " ax7.scatter(x[labels==0], y_corr[labels==0]/7.3,color='tab:blue')\n", + " ax7.scatter(x[labels==1], y_mgcorr/7.3,color='tab:orange')\n", + " ax7.scatter(x[labels==2], y_lgcorr/7.3,color='tab:green')\n", + " ax7.set_xscale('log')\n", + "\n", + " \n", + " else:\n", + " no_entry += 1\n", + "\n", + " ax6.hlines(np.mean(reldiff), 0, counter, color='k', zorder=10)\n", + " leg = ax5.legend()\n", + "\n", + " \n", + " fig1.show()\n", + " fig2.show()\n", + " fig3.show()\n", + " fig7.show()\n", + " print(f\"Pixels with signal not being injected in all gain stages: {no_entry//len(test_cells)}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tpix_range1 = [(200,210), (63,73)]\n", + "test_cells = [10,300]\n", + "roi = 'ROI1'\n", + "evaluate_fitting_ROI(cs_data[0], tpix_range1, test_cells, roi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tpix_range2 = [(35,38), (38,44)]\n", + "roi2 = 'ROI2'\n", + "evaluate_fitting_ROI(cs_data[0], tpix_range2, test_cells, roi2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def parallel_fit(yca, cor, variable):\n", + " \n", + " gain_marker = [0, 1, 2]\n", + " all_cores = multiprocessing.cpu_count()\n", + " num_workers = cor.shape[1] if cor.shape[1] < 50 else all_cores - 3\n", + " context= psh.context.ProcessContext(num_workers=num_workers)\n", + " \n", + " all_fres = {g: {var: context.alloc(shape=yca.shape[1:], \n", + " dtype=np.float32) for var in variable} for g in gain_marker}\n", + "\n", + " failures = context.alloc(shape=yca.shape[1:], dtype=np.int32) \n", + "\n", + " \n", + " def fit_row_parallel(worker_id, array_index, row):\n", + "\n", + " for cell in range(yca.shape[1]):\n", + " y = yca[:, cell, row]\n", + " x = x_eq[:y.shape[0]]\n", + " y_cor = cor[:, cell, row]\n", + " \n", + " try:\n", + " failures[cell, row] = 0\n", + " labels, gain_marker, bounds = label_gain_stage(x, y)\n", + "\n", + " if (np.any(bounds)):\n", + "\n", + " for g in gain_marker:\n", + " nonzero = np.where(y_cor[labels==g] > 0.)\n", + " x_sg = x[labels==g][nonzero][:bounds[g]]\n", + " if g == 0:\n", + " y_sg = y_cor[labels==g][nonzero][:bounds[g]]\n", + " else:\n", + " y_sg = y[labels==g][nonzero][:bounds[g]]\n", + "\n", + " if y_sg.shape[0] != 0: \n", + " parms = {'m': 1, 'b': np.min(y_sg)}\n", + " fitted = fit_data(lin_fun, x_sg, y_sg, parms)\n", + " yf = lin_fun(x_sg, fitted['m'], fitted['b'])\n", + "\n", + " all_fres[g]['m'][cell, row] = fitted['m']\n", + " all_fres[g]['b'][cell, row] = fitted['b']\n", + " all_fres[g]['error_m'][cell, row] = fitted['error_m']\n", + " all_fres[g]['red_chi2'][cell, row] = fitted['red_chi2'] \n", + " all_fres[g]['fit_dev'][cell, row] = np.median(np.abs((y_sg - yf) / y_sg))\n", + " all_fres[g]['no_entry'][cell, row] = 0 \n", + "\n", + " else:\n", + " all_fres[g]['no_entry'][cell, row] = 1\n", + "\n", + " except Exception as e:\n", + " failures[cell, row] = 1\n", + " \n", + " \n", + " context.map(fit_row_parallel, range(yca.shape[-1])) \n", + "\n", + " return all_fres, failures" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def execute_parallel_fit(data, corr_data):\n", + " gains = [0, 1, 2]\n", + " variable = ['m', 'b', 'error_m', 'red_chi2', 'fit_dev', 'no_entry']\n", + " fres = {g: {var: np.full(data.shape[1:], np.nan) for var in variable} for g in gains}\n", + "\n", + " failures_fit = np.full(data.shape[1:], np.nan)\n", + " \n", + " for col in range(data.shape[2]):\n", + " yca = data[..., col, :]\n", + " cor = corr_data[..., col, :]\n", + " frs = parallel_fit(yca, cor, variable)\n", + "\n", + " for gain in gains:\n", + "\n", + " fres[gain]['m'][:, col, :] = frs[0][gain]['m']\n", + " fres[gain]['b'][:, col, :]= frs[0][gain]['b']\n", + " fres[gain]['error_m'][:, col, :] = frs[0][gain]['error_m']\n", + " fres[gain]['red_chi2'][:, col, :] = frs[0][gain]['red_chi2'] \n", + " fres[gain]['fit_dev'][:, col, :] = frs[0][gain]['fit_dev']\n", + " fres[gain]['no_entry'][:, col, :] = frs[0][gain]['no_entry'] \n", + "\n", + " failures_fit[:, col, :] = frs[1]\n", + "\n", + " return fres, failures_fit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cor = cs_data[0] - offset" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true, + "tags": [] + }, + "outputs": [], + "source": [ + "start = datetime.now()\n", + "fres, failures_fit = execute_parallel_fit(cs_data[0], cor)\n", + "end = datetime.now() - start\n", + "print('Duration of fitting: ', end)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import copy\n", + "\n", + "fres_copy = copy.deepcopy(fres)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "failure_percentage = np.count_nonzero(failures_fit)/352/(128*512)*100\n", + "print('Percentage of fitting failure: {:.2f}'.format(failure_percentage))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "slopes = {}\n", + "intercepts = {}\n", + "red_chi2 = {}\n", + "for i, (gain, key) in enumerate(fres.items()):\n", + " slopes[i] = key['m']\n", + " intercepts[i] = key['b']\n", + " red_chi2[i] = key['red_chi2']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def calc_median(roi):\n", + " index = []\n", + " for idx in range(0, 128,roi[0]):\n", + " for j in range(0,(512//roi[1])):\n", + " index.append([[j*roi[1], j*roi[1]+roi[1]], [idx, idx+roi[0]]])\n", + " index = np.asarray(index)\n", + " \n", + " fshape = fres[0]['m'].shape\n", + " median_m = np.zeros((3,fshape[0], fshape[1], fshape[2]))\n", + " median_b = np.zeros((3,fshape[0], fshape[1], fshape[2]))\n", + "\n", + " for g in range(0,3):\n", + " means = np.zeros((fshape[0], fshape[1], fshape[2]))\n", + " means_b = np.zeros((fshape[0], fshape[1], fshape[2]))\n", + "\n", + " for r_ac in range(index.shape[0]):\n", + " idx = index[r_ac]\n", + " val = np.nanmedian(slopes[g][:, idx[1][0]:idx[1][1], \n", + " idx[0][0]:idx[0][1]], axis=(1,2))\n", + " means[:, idx[1][0]:idx[1][1], idx[0][0]:idx[0][1]]=np.repeat(val, roi[0]*roi[1]).reshape(fshape[0],roi[0],roi[1])\n", + "\n", + " val = np.nanmedian(intercepts[g][:, idx[1][0]:idx[1][1], \n", + " idx[0][0]:idx[0][1]], axis=(1,2))\n", + " means_b[:, idx[1][0]:idx[1][1], idx[0][0]:idx[0][1]]=np.repeat(val, roi[0]*roi[1]).reshape(fshape[0],roi[0],roi[1])\n", + " median_m[g, ...] = means\n", + " median_b[g, ...] = means_b\n", + " return median_m, median_b" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "roi = [4,8]\n", + "median_m, median_b = calc_median(roi)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example of high gain median values of cell 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = xana.heatmapPlot(median_m[0,1,...], cmap='jet', add_panels=False, aspect=1,\n", + " lut_label='Median slope', cb_loc='bottom', x_label='Rows', y_label='Columns')\n", + "fig = xana.heatmapPlot(median_b[0,1,...], cmap='jet', add_panels=False, aspect=1,\n", + " lut_label='Median intercept', cb_loc='bottom', x_label='Rows', y_label='Columns')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bad_pixels = {}\n", + "for g in range(0,3):\n", + " mask = np.zeros(slopes[g].shape, np.uint32)\n", + " error_m_thr = np.nanmean(fres[g]['error_m'])+5*np.nanstd(fres[g]['error_m'])\n", + " mask[(np.abs(fres[g]['error_m']) > error_m_thr) & (mask==0)] |= BadPixels.CI_EVAL_ERROR.value\n", + " mask[(np.abs(fres[g]['error_m']) == 0) & (mask==0)] |= BadPixels.CI_EVAL_ERROR.value\n", + " mask[(~np.isfinite(fres[g]['error_m'])) & (mask==0)] |= BadPixels.CI_EVAL_ERROR.value\n", + " mask[(fres[g]['red_chi2'] > chi2_lim) & (mask==0)] |= BadPixels.CI_EVAL_ERROR.value \n", + " mask[(~np.isfinite(fres[g]['red_chi2'])) & (mask==0)] |= BadPixels.CI_EVAL_ERROR.value \n", + " mask[(~np.isfinite(fres[g]['fit_dev'])) & (mask==0)] |= BadPixels.CI_LINEAR_DEVIATION.value\n", + " fit_dev_thr = np.nanmean(fres[g]['fit_dev'][np.isfinite(fres[g]['fit_dev'])])+5*np.nanstd(fres[g]['fit_dev'][np.isfinite(fres[g]['fit_dev'])]) \n", + " mask[(fres[g]['fit_dev'] > fit_dev_thr) & (mask==0)] |= BadPixels.CI_LINEAR_DEVIATION.value\n", + " mask[((slopes[g] <= 0.1) | (intercepts[g] <= 1.)) & (mask==0)] |= BadPixels.CI_EVAL_ERROR.value\n", + " mask[((~np.isfinite(slopes[g])) | (~np.isfinite(intercepts[g]))) & (mask==0)] |= BadPixels.CI_EVAL_ERROR.value \n", + "\n", + " bad_pixels[g] = mask" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fshape = fres[0]['m'].shape\n", + "BPmap = np.zeros(fshape, np.uint32) \n", + "for g in range(0,3):\n", + " BPmap[(bad_pixels[g] > 0) & (BPmap==0)] = bad_pixels[g][(bad_pixels[g] > 0) & (BPmap==0)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tshape = (fshape[0],fshape[1]*fshape[2])\n", + "\n", + "ratio_HM = slopes[0] / slopes [1]\n", + "ratio_ML = slopes[1] / slopes [2]\n", + "ratio_HM_med = median_m[0] / median_m[1]\n", + "ratio_ML_med = median_m[1] / median_m[2]\n", + "\n", + "rms_ratio_HM = mad(ratio_HM.reshape(tshape), axis=1, scale='normal', nan_policy='omit')\n", + "rms_ratio_ML = mad(ratio_ML.reshape(tshape), axis=1, scale='normal', nan_policy='omit')\n", + "rms_ratio_HM = np.repeat(rms_ratio_HM, fshape[1]*fshape[2]).reshape(fshape)\n", + "rms_ratio_ML = np.repeat(rms_ratio_ML, fshape[1]*fshape[2]).reshape(fshape)\n", + "\n", + "thr_mean_HM = [relgain_lim[0]*np.mean(ratio_HM_med[np.isfinite(ratio_HM_med)]), relgain_lim[1]*np.mean(ratio_HM_med[np.isfinite(ratio_HM_med)])]\n", + "thr_mean_ML = [relgain_lim[0]*np.mean(ratio_ML_med[np.isfinite(ratio_ML_med)]), relgain_lim[1]*np.mean(ratio_ML_med[np.isfinite(ratio_ML_med)])]\n", + "\n", + "bad_fits = (np.abs((ratio_HM - ratio_HM_med)/rms_ratio_HM) > sigma_dev_cut) |\\\n", + " (np.abs((ratio_ML - ratio_ML_med)/rms_ratio_ML) > sigma_dev_cut) |\\\n", + " (~np.isfinite(ratio_HM)) | (~np.isfinite(ratio_ML)) |\\\n", + " (ratio_HM < thr_mean_HM[0]) | (ratio_HM > thr_mean_HM[1]) |\\\n", + " (ratio_ML < thr_mean_ML[0]) | (ratio_ML > thr_mean_ML[1])\n", + "\n", + "BPmap[(bad_fits) & (BPmap==0)] |= BadPixels.CI_GAIN_OF_OF_THRESHOLD.value \n", + "bp_vals = [BadPixels.CI_GAIN_OF_OF_THRESHOLD.value, BadPixels.CI_LINEAR_DEVIATION.value, BadPixels.CI_EVAL_ERROR.value]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Bad pixel map" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nonzeros = []\n", + "for i in range(bad_fits.shape[0]):\n", + " nonzeros.append(np.count_nonzero(bad_fits[i]))\n", + "\n", + "worst_cell = np.where(nonzeros == np.max(nonzeros))[0][0] \n", + "print('Percentage of bad pixels on avg.: {:.2f}%'.format(np.mean(np.asarray(nonzeros))/(512*128)*100))\n", + "\n", + "fig=xana.heatmapPlot(bad_fits[1, ...], cmap='magma', add_panels=False, aspect=1,\n", + " lut_label='Bad pixels, cell: 1', cb_loc='bottom')\n", + "fig=xana.heatmapPlot(bad_fits[worst_cell, ...], cmap='magma', add_panels=False, aspect=1,\n", + " lut_label=f'Bad pixels (the worst cell), cell: {worst_cell}', cb_loc='bottom')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig=plt.figure(figsize=(6,4))\n", + "plt.plot(nonzeros, ls='--', marker='o')\n", + "plt.ylabel('# Bad pixels')\n", + "plt.xlabel('Memory cell Id')\n", + "plt.grid(axis='y', ls='dotted')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "total = np.count_nonzero(BPmap)\n", + "lbl = ['Gain off', 'Linear deviation', 'Evaluation error']\n", + "color = plt.cm.YlOrBr(np.linspace(0.2, 0.6, 3))\n", + "cut_flow = [np.count_nonzero(BPmap == 16), # CI_GAIN_OF_OF_THRESHOLD\n", + " np.count_nonzero((BPmap == 16) | (BPmap == 32)), #previous+CI_LINEAR_DEVIATION\n", + " np.count_nonzero(BPmap)] #previous + CI_EVAL_ERROR \n", + "fig, ax = plt.subplots(figsize=(8,4))\n", + "ax.yaxis.grid(zorder=0, ls='--')\n", + "lines = []\n", + "for i, bp_vl in enumerate(bp_vals):\n", + " ctr = np.count_nonzero(BPmap == bp_vl)\n", + " y = (1-ctr/np.prod(BPmap.shape))*100\n", + " y1 = (1-cut_flow[i]/np.prod(BPmap.shape))*100\n", + " lines += plt.bar(lbl[i], y, color='tab:blue', label='this cut', zorder=5, )\n", + " plt.bar(lbl[i], y1, color='tab:orange', alpha=1, zorder=5)\n", + "ax.set_xticklabels(lbl, rotation=45, ha='right', fontsize=12)\n", + "plt.ylim(y1-1, 100)\n", + "plt.ylabel('Percentage of good pixels (%)', fontsize=12)\n", + "labels = ['Single contribution', 'Cut flow']\n", + "plt.legend(labels)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def plot_cuts(a,sel,name, ax = None) :\n", + " \n", + " stats_cut = {}\n", + " stats = {}\n", + " if ax is None :\n", + " fig = plt.figure(figsize=(3,5))\n", + " ax = fig.add_subplot(111)\n", + " stats_cut[\"Mean\"] = np.nanmean(a[sel])\n", + " stats_cut[\"STD\"] = np.nanstd(a[sel])\n", + " stats_cut[\"Median\"] = np.nanmedian(a[sel])\n", + " stats_cut[\"Min\"] = np.nanmin(a[sel])\n", + " stats_cut[\"Max\"] = np.nanmax(a[sel])\n", + "\n", + " msk = (np.isfinite(a)) & (a > 0) & (a<14e3)\n", + " stats[\"Mean\"] = np.nanmean(a[msk])\n", + " stats[\"STD\"] = np.nanstd(a[msk])\n", + " stats[\"Median\"] = np.nanmedian(a[msk])\n", + " stats[\"Min\"] = np.nanmin(a[msk])\n", + " stats[\"Max\"] = np.nanmax(a[msk])\n", + "\n", + " m=np.nanmean(a[msk])\n", + " s=np.nanstd(a[msk])\n", + " bins = np.linspace(m-10*s,m+10*s,100)\n", + " ax.hist(a.ravel(), \n", + " bins = bins,\n", + " color='tab:red',\n", + " label='all fits',\n", + " alpha=0.7\n", + " )\n", + " ax.hist(a[sel].ravel(), \n", + " bins = bins,\n", + " color='tab:green',\n", + " label='good fits'\n", + " )\n", + " ax.grid(zorder=0)\n", + " ax.set_xlabel(name)\n", + " ax.set_yscale('log')\n", + "# ax.set_title(name)\n", + " ax.legend()\n", + " def statistic(stat, colour, shift):\n", + " textstr = \"\"\n", + " for key in stat:\n", + " try:\n", + " textstr += '{0: <6}: {1: 7.2f}\\n'.format(key, stat[key])\n", + " except:\n", + " pass\n", + " props = dict(boxstyle='round', alpha=0.5, color=colour,)\n", + " ax.text(0.05, 0.95-shift, textstr, transform=ax.transAxes, fontsize=10,\n", + " family='monospace', weight='book', stretch='expanded', \n", + " verticalalignment='top', bbox=props, zorder=2)\n", + " statistic(stats_cut, 'tab:green', 0)\n", + " statistic(stats, 'tab:red', 0.3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(15,15))\n", + "\n", + "plot_cuts(slopes[0],~(BPmap > 0), 'Slopes HG (ADU/DAC)', fig.add_subplot(321))\n", + "plot_cuts(slopes[1],~(BPmap > 0), 'Slopes MG (ADU/DAC)', fig.add_subplot(323))\n", + "plot_cuts(slopes[2],~(BPmap > 0), 'Slopes LG (ADU/DAC)', fig.add_subplot(325))\n", + "\n", + "plot_cuts(intercepts[0],~(BPmap > 0), 'Intercepts HG (ADU)', fig.add_subplot(322))\n", + "plot_cuts(intercepts[1],~(BPmap > 0), 'Intercepts MG (ADU)', fig.add_subplot(324))\n", + "plot_cuts(intercepts[2],~(BPmap > 0), 'Intercepts LG (ADU)', fig.add_subplot(326))\n", + "\n", + "fig = plt.figure(figsize=(15,15))\n", + "plot_cuts(ratio_HM,~(BPmap > 0), 'Ratio HG/MG', fig.add_subplot(321))\n", + "plot_cuts(ratio_ML,~(BPmap > 0), 'Ratio MG/LG', fig.add_subplot(322))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for g in range(0,3):\n", + " slopes[g][BPmap > 0] = median_m[g, ...][BPmap > 0]\n", + " intercepts[g][BPmap > 0] = median_b[g, ...][BPmap > 0]\n", + "ratio_HM = slopes[0] / slopes [1]\n", + "ratio_ML = slopes[1] / slopes [2]\n", + "\n", + "ratio_HM[np.isnan(ratio_HM) | np.isinf(ratio_HM)] = 0.\n", + "ratio_ML[np.isnan(ratio_ML) | np.isinf(ratio_ML)] = 0." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1,2, figsize=(10,5))\n", + "ax[0].hist(ratio_HM.flatten(), bins=100, range=thr_mean_HM, color='tab:green')\n", + "ax[1].hist(ratio_ML.flatten(), bins=100, range=thr_mean_ML, color='tab:green')\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from mpl_toolkits.axes_grid1 import AxesGrid\n", + "def preview_fitted_params(data, gain, cell_to_preview):\n", + " \n", + " g_label = ['High gain', 'Medium gain', 'Low gain']\n", + " \n", + " plot_text= {\"m\": 'Fitted slope (m)',\n", + " \"b\": 'Fitted intercept (b)',\n", + " \"error_m\": 'Slope error',\n", + " \"red_chi2\": 'Reduced Chi2',\n", + " \"fit_dev\": 'Relative difference fit vs. data',\n", + " \"no_entry\": 'No entry for fitting'\n", + " }\n", + " \n", + "\n", + " fig = plt.figure(figsize=(20,7))\n", + " plt.suptitle(f'{g_label[g]}', fontsize=16)\n", + "\n", + " grid = AxesGrid(fig, 111,\n", + " nrows_ncols=(3, 2),\n", + " axes_pad=(0.9, 0.55),\n", + " label_mode=\"0\",\n", + " share_all=True,\n", + " cbar_location=\"right\",\n", + " cbar_mode=\"each\",\n", + " cbar_size=\"7%\",\n", + " cbar_pad=\"2%\"\n", + " )\n", + " i = 0\n", + " \n", + " \n", + " for key, citem in data.items():\n", + " item = citem.copy()\n", + " item[~np.isfinite(item)] = 0\n", + " med = np.nanmedian(item)\n", + " bound = 0.1\n", + " maxcnt = 10\n", + " if med < 0:\n", + " bound = -bound\n", + "\n", + " while(np.count_nonzero((item < med-bound*med) | (item > med+bound*med))/item.size > 0.01): \n", + " bound *= 2\n", + " maxcnt -= 1\n", + " if maxcnt < 0:\n", + " break\n", + "\n", + " vmin = med-bound*med\n", + " vmax=med+bound*med\n", + "\n", + " if \"no_entry\" in key:\n", + " vmin = -1\n", + " vmax = 1\n", + " d = item[cell_to_preview, ...]\n", + " im = grid[i].imshow(d, interpolation=\"nearest\", origin='lower', cmap='rainbow',\n", + " vmin=vmin, vmax=vmax)\n", + " grid[i].set_title(plot_text[key], fontsize=14)\n", + " \n", + " cb = grid.cbar_axes[i].colorbar(im)\n", + " # axes coordinates are 0,0 is bottom left and 1,1 is upper right\n", + " x0, x1 = tpix_range1[0][0], tpix_range1[0][1]\n", + " y0, y1 = tpix_range1[1][0], tpix_range1[1][1]\n", + " p = patches.Rectangle((x0, y0), (x1 - x0), (y1 - y0), fill=False, color=\"red\")\n", + " grid[i].add_patch(p)\n", + " x0, x1 = tpix_range2[0][0], tpix_range2[0][1]\n", + " y0, y1 = tpix_range2[1][0], tpix_range2[1][1]\n", + " p = patches.Rectangle((x0, y0), (x1 - x0), (y1 - y0), fill=False, color=\"white\")\n", + " grid[i].add_patch(p)\n", + "\n", + " i += 1\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print('Results from cell: 1')\n", + "for g, (gain, data) in enumerate(fres.items()):\n", + " preview_fitted_params(data, g, 1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print('Results from cell: 200')\n", + "for g, (gain, data) in enumerate(fres.items()):\n", + "# print(g)\n", + " preview_fitted_params(data, g, 200)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "roi = [32,64]\n", + "median_m, median_b = calc_median(roi)\n", + "\n", + "for g in range(0,3):\n", + "# slopes[g][bad_fits] = median_m[g, ...][bad_fits]\n", + " slopes[g][BPmap > 0] = median_m[g, ...][BPmap > 0]\n", + "# intercepts[g][bad_fits] = median_b[g, ...][bad_fits]\n", + " intercepts[g][BPmap > 0] = median_b[g, ...][BPmap > 0]\n", + "ratio_HM = slopes[0] / slopes [1]\n", + "ratio_ML = slopes[1] / slopes [2]\n", + "\n", + "ratio_HM[np.isnan(ratio_HM) | np.isinf(ratio_HM)] = 0.\n", + "ratio_ML[np.isnan(ratio_ML) | np.isinf(ratio_ML)] = 0." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fres = copy.deepcopy(fres_copy) # this is needed to have raw fits without sanitization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if local_output:\n", + " for module in modules:\n", + " constants_file = '{}/CSconst_{}_M{}.h5'.format(out_folder, karabo_id, module)\n", + " with h5py.File(constants_file, 'w') as f:\n", + " for gain in fres:\n", + " for key, item in fres[gain].items():\n", + " f[f'/GainFits/{gain}/{key}/data'] = item\n", + " f['/BadPixels/data'] = BPmap\n", + " for g in range(0,3):\n", + " f[f'/SanitizedConsts/{g}/m/data'] = slopes[g]\n", + " f[f'/SanitizedConsts/{g}/b/data'] = intercepts[g]\n", + " f['/SanitizedConsts/Ratios/H-M/data'] = ratio_HM\n", + " f['/SanitizedConsts/Ratios/M-L/data'] = ratio_ML\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "## INJECTION OF CS CONSTS TO DB DOES NOT WORK YET!\n", + "\n", + "# md = None\n", + "\n", + "# # set the operating condition\n", + "# condition = Conditions.Dark.AGIPD(memory_cells=mem_cells, \n", + "# bias_voltage=bias_voltage,\n", + "# acquisition_rate=acq_rate, \n", + "# gain_setting=gain_setting,\n", + "# integration_time=integration_time)\n", + "\n", + "# db_modules = get_pdu_from_db(karabo_id, karabo_da, Constants.AGIPD.SlopesPC(),\n", + "# condition, cal_db_interface,\n", + "# snapshot_at=creation_time)\n", + "\n", + "# for pdu, (qm, r) in zip(db_modules, fres.items()):\n", + "# for const in [\"SlopesCS\", \"BadPixelsCS\"]:\n", + "\n", + "# dbconst = getattr(Constants.AGIPD, const)()\n", + "\n", + "# if const == \"SlopesPC\":\n", + "# dbconst.data = slope_dict_to_arr(r)\n", + "# else:\n", + "# dbconst.data = bad_pixels[qm]\n", + "\n", + "# if db_output:\n", + "# md = send_to_db(pdu, karabo_id, dbconst, condition,\n", + "# file_loc, report, cal_db_interface,\n", + "# creation_time=creation_time,\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\"• 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": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print('Ratios from cell: 1')\n", + "cell = 1\n", + "vmin = 0.9*np.min(ratio_HM[cell, ...][ratio_HM[cell, ...]>0.])\n", + "vmax = 1.3*np.mean(ratio_HM[cell, ...][ratio_HM[cell, ...]>0.])\n", + "fig=xana.heatmapPlot(ratio_HM[cell, ...], cmap='jet', add_panels=False, aspect=1,\n", + " lut_label='Cell: {} slope ratio HG / MG'.format(cell), cb_loc='bottom', vmin=vmin, vmax=vmax)\n", + "\n", + "vmin = 0.9*np.min(ratio_ML[cell, ...][ratio_ML[cell, ...]>0.])\n", + "vmax = 1.3*np.mean(ratio_ML[cell, ...][ratio_ML[cell, ...]>0.])\n", + "fig=xana.heatmapPlot(ratio_ML[cell, ...], cmap='jet', add_panels=False, aspect=1,\n", + " lut_label='Cell: {} slope ratio MG / LG'.format(cell), cb_loc='bottom', vmin=vmin, vmax=vmax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print('Ratios from cell: 200')\n", + "cell = 200\n", + "vmin = 0.9*np.min(ratio_HM[cell, ...][ratio_HM[cell, ...]>0.])\n", + "vmax = 1.3*np.mean(ratio_HM[cell, ...][ratio_HM[cell, ...]>0.])\n", + "fig=xana.heatmapPlot(ratio_HM[cell, ...], cmap='jet', add_panels=False, aspect=1,\n", + " lut_label='Cell: {} slope ratio HG / MG'.format(cell), cb_loc='bottom', vmin=vmin, vmax=vmax)\n", + "\n", + "vmin = 0.9*np.min(ratio_ML[cell, ...][ratio_ML[cell, ...]>0.])\n", + "vmax = 1.3*np.mean(ratio_ML[cell, ...][ratio_ML[cell, ...]>0.])\n", + "fig=xana.heatmapPlot(ratio_ML[cell, ...], cmap='jet', add_panels=False, aspect=1,\n", + " lut_label='Cell: {} slope ratio MG / LG'.format(cell), cb_loc='bottom', vmin=vmin, vmax=vmax)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Mean value of gain ratios" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "HGMG_mean = np.mean(ratio_HM, axis=0)\n", + "MGLG_mean = np.mean(ratio_ML, axis=0)\n", + "MGLG_mean[~np.isfinite(MGLG_mean)] = np.median(MGLG_mean)\n", + "\n", + "vmin = 0.9*np.mean(HGMG_mean[HGMG_mean > 0.])\n", + "vmax = 1.2*np.mean(HGMG_mean[HGMG_mean > 0.])\n", + "# print(vmin, vmax)\n", + "fig=xana.heatmapPlot(HGMG_mean, cmap='jet', add_panels=False, aspect=1,\n", + " lut_label='Mean slope ratio HG / MG', cb_loc='bottom', vmin=vmin, vmax=vmax)\n", + "\n", + "vmin = 0.9*np.mean(MGLG_mean[MGLG_mean > 0.])\n", + "vmax = 1.2*np.mean(MGLG_mean[MGLG_mean > 0.])\n", + "fig=xana.heatmapPlot(MGLG_mean, cmap='jet', add_panels=False, aspect=1,\n", + " lut_label='Mean slope ratio MG / LG', cb_loc='bottom', vmin=vmin, vmax=vmax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(3, 2, figsize=(16, 8))\n", + "fig.subplots_adjust(wspace=0.1, hspace=0.5)\n", + "lbl = ['HG', 'MG', 'LG']\n", + "for i in range(3):\n", + " hrange = [0.1 * np.median(slopes[i]), 2.5 * np.median(slopes[i])]\n", + " r_hist = np.zeros((slopes[i].shape[0], 100))\n", + " for c in range(slopes[i].shape[0]):\n", + " d = slopes[i][c, ...]\n", + "\n", + " h, e = np.histogram(d.flatten(), bins=100, range=hrange)\n", + " r_hist[c, :] = h\n", + "\n", + " im = ax[i, 0].imshow(r_hist[:, :].T[::-1, :], interpolation=\"nearest\",\n", + " aspect=\"auto\", norm=LogNorm(vmin=1, vmax=np.max(r_hist)),\n", + " extent=[0, slopes[i].shape[0], hrange[0], hrange[1]])\n", + " ax[i, 0].set_xlabel(\"Memory cell\", fontsize=12)\n", + " ax[i, 0].set_ylabel('Slope {}'.format(lbl[i]), fontsize=12)\n", + " cb = fig.colorbar(im, ax=ax[i, 0], pad=0.01)\n", + " \n", + " hrange = [0.5 * np.median(intercepts[i]), 1.8 * np.median(intercepts[i])]\n", + " r_hist = np.zeros((intercepts[i].shape[0], 100))\n", + " for c in range(intercepts[i].shape[0]):\n", + " d = intercepts[i][c, ...]\n", + "\n", + " h, e = np.histogram(d.flatten(), bins=100, range=hrange)\n", + " r_hist[c, :] = h\n", + "\n", + " im = ax[i, 1].imshow(r_hist[:, :].T[::-1, :], interpolation=\"nearest\",\n", + " aspect=\"auto\", norm=LogNorm(vmin=1, vmax=np.max(r_hist)),\n", + " extent=[0, intercepts[i].shape[0], hrange[0], hrange[1]])\n", + " ax[i, 1].set_xlabel(\"Memory cell\", fontsize=12)\n", + " ax[i, 1].set_ylabel('Intercept {}'.format(lbl[i]), fontsize=12)\n", + " cb = fig.colorbar(im, ax=ax[i, 1], pad=0.01, label='Frequency')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "ratios = [ratio_HM, ratio_ML]\n", + "fig, ax = plt.subplots(1, 2, figsize=(16, 4))\n", + "lbl = ['High-Medium', 'Medium-Low']\n", + "for i, ratio in enumerate(ratios):\n", + " hrange = [0.5 * np.median(ratio), 1.5 * np.median(ratio)]\n", + "# print(hrange)\n", + " r_hist = np.zeros((ratio.shape[0], 100))\n", + " for c in range(ratio.shape[0]):\n", + " d = ratio[c, ...]\n", + "\n", + " h, e = np.histogram(d.flatten(), bins=100, range=hrange)\n", + " r_hist[c, :] = h\n", + "\n", + " im = ax[i].imshow(r_hist[:, :].T[::-1, :], interpolation=\"nearest\",\n", + " aspect=\"auto\", norm=LogNorm(vmin=1, vmax=np.max(r_hist)),\n", + " extent=[0, ratio_HM.shape[0], hrange[0], hrange[1]])\n", + " ax[i].set_xlabel(\"Memory cell\", fontsize=12)\n", + " ax[i].set_ylabel('Ratio {}'.format(lbl[i]), fontsize=12)\n", + " cb = fig.colorbar(im, ax=ax[i], pad=0.01, label='Frequency')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.11" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/notebooks/AGIPD/CS_parallelMerging_NBC.ipynb b/notebooks/AGIPD/CS_parallelMerging_NBC.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..63f6073a34f49270df1db17b267a42c6cb66d4e3 --- /dev/null +++ b/notebooks/AGIPD/CS_parallelMerging_NBC.ipynb @@ -0,0 +1,604 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Merging of Current Source Data\n", + "Author: Detector Group, Version 1.0\n", + "\n", + "This notebook is used to create for each AGIPD module a single .h5 file containing current source injected data. Current source (CS) data for each module consists of four runs as at a time the current source can injected charge only in every 4th column. The merged .h5 file is then used by notebook \"CS_Characterization_unequalClockStep_NBC\" to create CS constants." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# cluster_profile = \"noDB\" # The ipcluster profile to use\n", + "in_folder='/gpfs/exfel/exp/SPB/202231/p900295/raw/' # path to input data, required\n", + "out_folder = \"/gpfs/exfel/exp/SPB/202230/p900250/scratch/CSmergedFiles\" # path to output to, required\n", + "runs1 = [19, 20, 21, 22] # runs to use, required, range allowed ITESTC 65\n", + "runs2 = [14, 15, 16, 17] # ITESTC 100\n", + "runs3 = [27, 28, 29, 30] # ITESTC 150\n", + "modules = [15] # modules to work on, required, range allowed\n", + "\n", + "karabo_da = [\"all\"]\n", + "karabo_da_control = \"AGIPD1MCTRL00\" # karabo DA for control infromation\n", + "karabo_id_control = \"SPB_IRU_AGIPD1M1\" # karabo-id for the control device e.g. \"MID_EXP_AGIPD1M1\", or \"SPB_IRU_AGIPD1M1\"\n", + "karabo_id = \"SPB_DET_AGIPD1M-1\"\n", + "ctrl_source_template = '{}/MDL/FPGA_COMP' # path to control information\n", + "instrument_source_template = '{}/DET/{}:xtdf' # path in the HDF5 file to images\n", + "receiver_template = \"{}CH0\" # inset for receiver devices\n", + "\n", + "use_dir_creation_date = True\n", + "delta_time = 0 # offset to the creation time (e.g. useful in case we want to force the system to use diff. dark constants)\n", + "cal_db_interface = \"tcp://max-exfl016:8019\" # the database interface to use\n", + "\n", + "\n", + "bias_voltage = 300 # detector bias voltage\n", + "mem_cells = 352 # number of memory cells used, use -1 to auto-derive\n", + "acq_rate = 4.5 # the detector acquisition rate, use -1. to try to auto-determine\n", + "gain_setting = 0 # gain setting can have value 0 or 1, -1 for auto-detection\n", + "integration_time = 12 # integration time, negative values for auto-detection. \n", + "\n", + "steps = [1, 10, 50]\n", + "increments = [300, 500, 100]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os, psutil\n", + "import warnings\n", + "from datetime import datetime, timedelta\n", + "from functools import partial\n", + "\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "import dateutil.parser\n", + "import h5py\n", + "import matplotlib\n", + "import numpy as np\n", + "\n", + "from scipy.stats import median_abs_deviation as mad\n", + "from collections import OrderedDict\n", + "import xarray\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import gridspec\n", + "\n", + "import XFELDetAna.xfelpyanatools as xana\n", + "from extra_data import RunDirectory, components\n", + "import pasha as psh\n", + "import multiprocessing\n", + "import itertools\n", + "from concurrent.futures import ProcessPoolExecutor\n", + "\n", + "from cal_tools.agipdlib import AgipdCtrl\n", + "from cal_tools.enums import BadPixels\n", + "from cal_tools.plotting import plot_badpix_3d, show_overview\n", + "from cal_tools.tools import (\n", + " gain_map_files,\n", + " get_constant_from_db_and_time,\n", + " get_dir_creation_date,\n", + " get_notebook_name,\n", + " get_pdu_from_db,\n", + " get_report,\n", + " module_index_to_qm,\n", + " parse_runs,\n", + " send_to_db,\n", + ")\n", + "from iCalibrationDB import Conditions, ConstantMetaData, Constants, Detectors, Versions\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cells = mem_cells\n", + "path_temp = in_folder+\"/r{:04d}/\"\n", + "image_name_temp = 'RAW-R{:04d}-AGIPD{:02d}-S{:05d}.h5'\n", + "print(\"Parameters are:\")\n", + "if mem_cells < 0:\n", + " print(\"Memory cells: auto-detection on\")\n", + "else:\n", + " print(\"Memory cells set by user: {}\".format(mem_cells))\n", + "print(\"Runs: {} / {} / {}\".format(runs1, runs2, runs3))\n", + "print(\"Modules: {}\".format(modules))\n", + "\n", + "instrument = karabo_id.split(\"_\")[0]\n", + "\n", + "\n", + "if instrument == \"HED\":\n", + " nmods = 8\n", + "else:\n", + " nmods = 16\n", + "\n", + "print(f\"Detector in use is {karabo_id}\")\n", + "\n", + "\n", + "if karabo_da == [\"all\"]:\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]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for module in modules:\n", + " first_run = runs2[0]\n", + " channel = module\n", + " ctrl_src = ctrl_source_template.format(karabo_id_control)\n", + " instrument_src = instrument_source_template.format(karabo_id, receiver_template).format(channel)\n", + "\n", + " agipd_cond = AgipdCtrl(\n", + " run_dc=RunDirectory(f'{in_folder}/r{first_run:04d}'),\n", + " image_src=instrument_src,\n", + " ctrl_src=ctrl_src,\n", + " raise_error=False, # to be able to process very old data without gain_setting value\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define creation time\n", + "creation_time = None\n", + "if use_dir_creation_date:\n", + " creation_time = get_dir_creation_date(in_folder, first_run)\n", + " creation_time = creation_time + timedelta(hours=delta_time)\n", + "\n", + "# Read AGIPD parameter conditions.\n", + "if acq_rate < 0.:\n", + " acq_rate = agipd_cond.get_acq_rate()\n", + "if mem_cells < 0:\n", + " mem_cells = agipd_cond.get_num_cells()\n", + "# TODO: look for alternative for passing creation_time\n", + "if gain_setting < 0:\n", + " gain_setting = agipd_cond.get_gain_setting(creation_time)\n", + "if bias_voltage < 0:\n", + " bias_voltage = agipd_cond.get_bias_voltage(karabo_id_control)\n", + "if integration_time < 0:\n", + " integration_time = agipd_cond.get_integration_time()\n", + "\n", + "print(f\"Acquisition rate: {acq_rate} MHz\")\n", + "print(f\"Memory cells: {mem_cells}\")\n", + "print(f\"Gain setting: {gain_setting}\")\n", + "print(f\"Integration time: {integration_time}\")\n", + "print(f\"Using {creation_time} as creation time of constant.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def count_min_bursts(in_folder, channel, run):\n", + " \n", + " bursts = []\n", + " print(run)\n", + " run_folder_path = in_folder.format(run)\n", + " r = RunDirectory(run_folder_path)\n", + " instrument_source = karabo_id+'/DET/{}CH0:xtdf'.format(channel)\n", + " print(instrument_source)\n", + " c = r.get_array(instrument_source, 'image.length')\n", + " total_frame = np.count_nonzero(c)\n", + " cells = r.detector_info(instrument_source)['frames_per_train']\n", + " bursts = total_frame // cells\n", + "\n", + " del c\n", + "\n", + " return bursts\n", + "\n", + "trains = []\n", + "for module in modules:\n", + " partial_check = partial(count_min_bursts, in_folder+\"/r{:04d}/\", module)\n", + "\n", + " from concurrent.futures import ProcessPoolExecutor\n", + "\n", + " with ProcessPoolExecutor(max_workers=4) as pool:\n", + " bursts1 = list(pool.map(partial_check, runs1))\n", + " bursts2 = list(pool.map(partial_check, runs2))\n", + " bursts3 = list(pool.map(partial_check, runs3)) \n", + " trains.append(np.min(bursts1 + bursts2 + bursts3))\n", + "trains = np.min(np.asarray(trains))\n", + "print('Number of trains: ',trains)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# check_ASIC looks which set of runs to use for a given ASIC\n", + "def check_ASIC(in_folder, runs, trains, channel, tresh):\n", + " \n", + " ASIC_code = []\n", + " \n", + " run_folder_path = in_folder.format(runs[0])\n", + " r = RunDirectory(run_folder_path)\n", + "\n", + " instrument_source = karabo_id+'/DET/{}CH0:xtdf'.format(channel) \n", + "\n", + " d = r.get_array(instrument_source, 'image.data')[:trains*cells, 0, :, :]\n", + " d = d.values.reshape((d.shape[0] // cells, cells, 1, 512, 128))\n", + " d = np.moveaxis(d, 4, 3).astype(np.float32)\n", + "\n", + " \n", + " for col in range(0, 512, 64):\n", + " rise_val = np.mean(d[-10:, 3, 0, :64, 3+col:64+col:4]) -\\\n", + " np.mean(d[800:810, 3, 0, :64, 3+col:64+col:4])\n", + " print(rise_val)\n", + " if (rise_val < tresh) & (rise_val > -30.):\n", + " ASIC_code.append(2)\n", + " else:\n", + " ASIC_code.append(1)\n", + " \n", + " for col in range(0, 512, 64): \n", + " rise_val = np.mean(d[-10:, 3, 0, 64:, col:64+col:4]) -\\\n", + " np.mean(d[800:810, 3, 0, 64:, col:64+col:4])\n", + " print(rise_val)\n", + " if (rise_val < tresh) & (rise_val > -30.):\n", + " ASIC_code.append(2)\n", + " else:\n", + " ASIC_code.append(1)\n", + " \n", + " del d\n", + " return ASIC_code\n", + "\n", + "\n", + "# Specify the runs to be used and ASIC indexing within data array in OrderedDict\n", + "def defineInitParams(m_num):\n", + " \n", + " channel = modules[m_num]\n", + " runs_dict = OrderedDict()\n", + " print(trains)\n", + " \n", + " instrument_source = karabo_id+'/DET/{}CH0:xtdf'.format(channel)\n", + " print(instrument_source)\n", + " \n", + " run_lbl = [0, runs1, runs2, runs3]\n", + " \n", + " for idx, a_c in enumerate(ASIC_code_list[m_num]):\n", + " \n", + " print('ASIC code:',a_c)\n", + " \n", + " if idx < 8:\n", + " index = [[idx*64, idx*64 + 64], [0, 64]]\n", + " else:\n", + " counter = 512\n", + " index = [[idx*64-counter, idx*64 + 64 - counter], [64, 128]]\n", + " counter -= 64\n", + " \n", + " run_cllection = [RunDirectory(f'{in_folder}/r{r_num:04d}/') for r_num in run_lbl[a_c]]\n", + " print('Used runs: ',run_lbl[a_c])\n", + " sel = [run_cllect.select(instrument_source, require_all=True) for run_cllect in run_cllection]\n", + " runs_dict[idx] = {'runs': sel,\n", + " 'module': channel,\n", + " 'index': index\n", + " }\n", + "\n", + " return runs_dict" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# %%time\n", + "# partial_check = partial(check_ASIC, in_folder+\"/r{:04d}/\", runs2[0], trains)\n", + "\n", + "# from concurrent.futures import ProcessPoolExecutor\n", + "\n", + "# with ProcessPoolExecutor(max_workers=2) as pool:\n", + "# ASIC_code = list(pool.map(partial_check, modules))\n", + "ASIC_code_list = []\n", + "for module in modules:\n", + " start = datetime.now()\n", + "\n", + " ASIC_code = check_ASIC(in_folder+\"/r{:04d}/\", runs1, trains, module, 500.)\n", + " print(ASIC_code)\n", + "\n", + " if 2 in ASIC_code:\n", + " print('\\n')\n", + " ASIC_code = np.asarray(ASIC_code)\n", + " ASIC_code2 = check_ASIC(in_folder+\"/r{:04d}/\", runs2, trains, module, 100.)\n", + " ASIC_code2 = np.array(ASIC_code2)\n", + " ASIC_code2[np.where(ASIC_code == 1)[0]] = 1\n", + " if 2 in ASIC_code2:\n", + "# ASIC_code2 = np.asarray(ASIC_code2)\n", + " substitute = np.where(ASIC_code2 == 2)[0]\n", + " ASIC_code[substitute] = 3\n", + " ASIC_code_list.append(ASIC_code)\n", + " end = datetime.now() - start\n", + " print('\\nDuration of check: ', end)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ASIC_code_list" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# with multiprocessing.Pool(processes=len(modules)) as pool:\n", + "# runs_dict_list = pool.map(defineInitParams, range(len(modules)))\n", + "runs_dict_list = []\n", + "for module in range(len(modules)):\n", + " runs_dict_list.append(defineInitParams(module))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def process_module(channel, runs_dict, cell_id):\n", + "\n", + " instrument_source = karabo_id+'/DET/{}CH0:xtdf'.format(runs_dict[channel][0]['module'])\n", + " print(instrument_source)\n", + "\n", + " context= psh.context.ProcessContext(num_workers=4)\n", + " cs_data = {'analog': context.alloc(shape=(trains, 352, 128, 512), dtype=np.float32),\n", + "# 'cellId': context.alloc(shape=(trains, 176), dtype=np.uint32)\n", + " # do we have to save also digital data?\n", + "# 'digital': context.alloc(shape=(trains, 36, 128, 512), dtype=np.float32),\n", + " }\n", + " \n", + " \n", + " def fill_ASIC(worker_id, array_index, asic):\n", + " run_asic = runs_dict[channel][asic]['runs']\n", + " module = runs_dict[channel][asic]['module']\n", + " index = runs_dict[channel][asic]['index']\n", + " print(\"Module: {}, ASIC: {}, Index: {}\".format(module, asic, index))\n", + " print('Cells: {} - {}'.format(cell_id[0], cell_id[1]))\n", + "\n", + " for counter, r in enumerate(run_asic):\n", + " d = r.get_array(instrument_source, 'image.data')[:trains*352, 0, index[0][0]:index[0][1], index[1][0]:index[1][1]]\n", + " d = d.values.reshape((d.shape[0] // mem_cells, mem_cells, 1, 512//8, 128//2))\n", + " d = np.moveaxis(d, 4, 3)\n", + " d = xarray.DataArray(d, coords={'cell': range(0, 352)}, dims=[\"train\", \"cell\", \"output\", \"x\", \"y\"])\n", + " d = d.sel(cell=range(cell_id[0], cell_id[1])) # select chunk of mem cells\n", + "\n", + "# cs_data['cellId'] = np.tile(d.cell.values, trains).reshape(trains, 44)\n", + " if asic < 8:\n", + " for i in range(3-counter, 64, 4):\n", + " cs_data['analog'][..., index[1][0]:index[1][1], index[0][0]+i] = d[:, :, 0, :, i]\n", + " else:\n", + " for i in range(counter, 64, 4):\n", + " cs_data['analog'][..., index[1][0]:index[1][1], index[0][0]+i] = d[:, :, 0, :, i]\n", + " del d\n", + " \n", + " context.map(fill_ASIC, range(0,16))\n", + " \n", + " return modules[channel], cs_data['analog']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "\n", + "run_lbl = [0, runs1, runs2, runs3]\n", + "\n", + "for module in range(len(modules)):\n", + " channel = modules[module]\n", + " values, counts = np.unique(ASIC_code_list[module], return_counts=True)\n", + " apendix = []\n", + " for val in values:\n", + " apendix.append('r{}-{}'.format(run_lbl[val][0], run_lbl[val][-1]))\n", + "# print(apendix)\n", + " merged_file = \"{}/agipd_CH{:02d}_cs_\".format(out_folder, channel)+\"_\".join(apendix)+'.h5'\n", + " print('Saving to:',merged_file)\n", + " \n", + " with h5py.File(merged_file, \"w\") as f:\n", + " dset = f.require_dataset(\"{}/Analog/data\".format(channel), (trains, mem_cells, 128, 512), dtype='float32')\n", + "# cellset = f.require_dataset(\"{}/CellId/data\".format(channel), (trains, mem_cells), dtype='uint32')\n", + " for cell in range(0,352,352):\n", + " cell_id = [cell, cell+352]\n", + " cs_data = process_module(module, runs_dict_list, cell_id) #cs_data[0] - channel, cs_data[1] - analog, cs_data[2] digital\n", + " dset[:, cell_id[0]:cell_id[1], ...] = cs_data[1]\n", + "# cellset[:, cell_id[0]:cell_id[1]] = cs_data[2]\n", + " f.flush()\n", + " del cs_data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Checking the resulting data after merging" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cs_data = []\n", + "with h5py.File(merged_file, 'r') as f:\n", + " cs_data.append(np.array(f[f'/{modules[0]}/Analog/data']))\n", + "\n", + "trains = cs_data[0].shape[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mc = 300\n", + "burst_of_interest = [10, 300, -50]\n", + "im_ratio = cs_data[0][0, 0, ...].shape[1] / cs_data[0][0, 0, ...].shape[0]\n", + "vmin = 0.9*np.mean(cs_data[0][mc, ...])\n", + "vmax = 14000\n", + "\n", + "fig, axs = plt.subplots(len(burst_of_interest), 1, figsize=(10,10))\n", + "fig.subplots_adjust(hspace=0.5)\n", + "yticks_major=np.arange(0,129,64)\n", + "xticks_major=np.arange(0,512,64)\n", + "for i, ax in enumerate(axs):\n", + " im = ax.imshow(cs_data[0][burst_of_interest[i], mc, ...], \n", + " cmap='jet', vmin=vmin, vmax=vmax)\n", + " cbar = fig.colorbar(im, ax=ax,fraction=0.0355*im_ratio, \n", + " pad=0.15, orientation='horizontal')\n", + " cbar.set_label('CS signal (ADU)')\n", + " b = burst_of_interest[i]\n", + " if burst_of_interest[i] < 0:\n", + " b = trains + b\n", + " ax.set_title(\"Scan point: {}, Cell: {}\".format(b, mc))\n", + " ax.set_xticks(xticks_major)\n", + " ax.xaxis.grid(True, which='major', linewidth=0.5, color='grey')\n", + " ax.set_yticks(yticks_major)\n", + " ax.yaxis.grid(True, which='major', linewidth=0.5, color='grey')\n", + " ax.set_xlabel('Columns')\n", + " ax.set_ylabel('Rows')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a = np.arange(0,increments[0],steps[0])\n", + "b = np.arange(a[-1]+1,a[-1]+1+increments[1]*steps[1], steps[1])\n", + "c = np.arange(b[-1]+1,b[-1]+1+increments[2]*steps[2], steps[2])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The following plots show integrated signal over the scan for a singel pixel in each ASIC for a memory cell # 4." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "row = 81\n", + "diff = np.sum(increments)-trains\n", + "mc = 300\n", + "\n", + "fig, axes = plt.subplots(1, 8, figsize=(14, 2), sharey=True, constrained_layout=True) \n", + "for asic, col in enumerate(range(35,485,64), start=8):\n", + " axes[asic-8].plot(a, cs_data[0][:increments[0],mc,row,col], ls=\"None\", marker='.', label='step size: 1 clk')\n", + " axes[asic-8].plot(b, cs_data[0][increments[0]:np.sum(increments[:2]),mc,row,col], ls=\"None\", marker='.', label='step size: 25 clk')\n", + " axes[asic-8].plot(c[:-diff], cs_data[0][np.sum(increments[:2]):,mc,row,col], ls=\"None\", marker='.', label='step size: 100 clk')\n", + " axes[asic-8].set_xlabel(asic)\n", + " axes[asic-8].set_xscale('log')\n", + "for ax in axes:\n", + " ax.set_xticks([])\n", + " ax.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "row = 31\n", + "diff = np.sum(increments)-trains\n", + "mc = 300\n", + "\n", + "fig, axes = plt.subplots(1, 8, figsize=(14, 2), sharey=True, constrained_layout=True) \n", + "for asic, col in enumerate(range(35,485,64)):\n", + " axes[asic].plot(a, cs_data[0][:increments[0],mc,row,col], ls=\"None\", marker='.', label='step size: 1 clk')\n", + " axes[asic].plot(b, cs_data[0][increments[0]:np.sum(increments[:2]),mc,row,col], ls=\"None\", marker='.', label='step size: 25 clk')\n", + " axes[asic].plot(c[:-diff], cs_data[0][np.sum(increments[:2]):,mc,row,col], ls=\"None\", marker='.', label='step size: 100 clk')\n", + " axes[asic].set_xlabel(asic)\n", + " axes[asic].set_xscale('log')\n", + "for ax in axes:\n", + " ax.set_xticks([])\n", + " ax.grid()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.11" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/src/xfel_calibrate/notebooks.py b/src/xfel_calibrate/notebooks.py index aaee10cacb262c48448b60be3bcbb5c60ff4761b..2a45d7264f6eb8d377fd163ec6c606a9b4af4462 100644 --- a/src/xfel_calibrate/notebooks.py +++ b/src/xfel_calibrate/notebooks.py @@ -18,6 +18,18 @@ notebooks = { "default concurrency": 16, "cluster cores": 32}, }, + "CS_MERGE": { + "notebook": "notebooks/AGIPD/CS_parallelMerging_NBC.ipynb", + "concurrency": {"parameter": "modules", + "default concurrency": 16, + "cluster cores": 1}, + }, + "CS": { + "notebook": "notebooks/AGIPD/CS_Characterization_unequalClockStep_NBC.ipynb", + "concurrency": {"parameter": "modules", + "default concurrency": 16, + "cluster cores": 1}, + }, "FF": { "notebook": "notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb",