{
"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/202331/p900376/raw/\" # path to input data, required\n",
"out_folder = \"/gpfs/exfel/exp/SPB/202331/p900376/scratch/CSmergedFiles/\" # path to output to, required\n",
"metadata_folder = \"\"\n",
"\n",
"first_run = 269 # first taken run, it has to be a run with the smallest ITESTC, otherwise define runs manually\n",
"runs1 = [-1] # list of runs to use, range allowed ITESTC 65, use -1 for auto-completion\n",
"runs2 = [-1] # list of runs to use, ITESTC 80, use -1 for auto-completion\n",
"runs3 = [-1] # list of runs to use, ITESTC 120, use -1 for auto-completion\n",
"runs4 = [-1] # list of runs to use, ITESTC 170, use -1 for auto-completion\n",
"itestc_order = \"descending\" # order in which ITESTC values were configured, e.g. first run taken with the highest ITESTC value --> \"descending\"\n",
"n_itestc = 5 # number of different ITESTC configurations taken. IMPORTANT, if you use ITESTC configurations which do not have consecutive run numbers, do not use runs auto-completion! Provide run numbers manually!\n",
"\n",
"modules = [-1] # modules to work on, required, range allowed\n",
"karabo_da = [\"all\"]\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",
"creation_time = \"\" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC e.g. \"2022-06-28 13:00:00\"\n",
"creation_date_offset = \"00:00:00\" # add an offset to creation date, e.g. to get different constants\n",
"cal_db_interface = \"tcp://max-exfl-cal001:8015#8045\" # the database interface to use\n",
"\n",
"bias_voltage = -1 # detector bias voltage\n",
"mem_cells = -1 # number of memory cells used, use -1 to use value stored in slow data.\n",
"acq_rate = -1. # the detector acquisition rate, use -1. to use value stored in slow data.\n",
"gain_setting = -1 # gain setting can have value 0 or 1, use -1 to use value stored in slow data.\n",
"integration_time = -1 # integration time, use -1 to use value stored in slow data.\n",
"\n",
"steps = [1, 10, 75] # spacing between integration time steps for each loop\n",
"increments = [300, 400, 200] # number of steps within a loop"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Determine the run numbers\n",
"if (runs1[0] == -1) & (itestc_order == \"ascending\"):\n",
" runs1 = range(first_run, first_run+4*4, n_itestc)\n",
" runs2 = range(first_run+1, first_run+1+4*4, n_itestc)\n",
" runs3 = range(first_run+2, first_run+2+4*4, n_itestc)\n",
" runs4 = range(first_run+3, first_run+3+4*4, n_itestc)\n",
" \n",
"if (runs1[0] == -1) & (itestc_order == \"descending\"):\n",
" runs1 = range(first_run, first_run+4*4, n_itestc)\n",
" runs2 = range(first_run-1, first_run-1+4*4, n_itestc)\n",
" runs3 = range(first_run-2, first_run-2+4*4, n_itestc)\n",
" runs4 = range(first_run-3, first_run-3+4*4, n_itestc)\n",
" \n",
"print('Runs to be evaluated: ', list(runs1), list(runs2), list(runs3), list(runs4))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"from datetime import datetime, timedelta\n",
"from functools import partial\n",
"from dateutil import parser\n",
"\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"import h5py\n",
"import matplotlib\n",
"import numpy as np\n",
"\n",
"from collections import OrderedDict\n",
"import xarray\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from extra_data import RunDirectory\n",
"import pasha as psh\n",
"from concurrent.futures import ProcessPoolExecutor\n",
"\n",
"from cal_tools.agipdlib import AgipdCtrl\n",
"from cal_tools.tools import calcat_creation_time\n",
"from cal_tools.step_timing import StepTimer\n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raw_dc = RunDirectory(in_folder+f'/r{first_run:04d}')\n",
"if karabo_da == [\"all\"]:\n",
" if modules[0] == -1:\n",
" strings = [s.split(\"/\")[-1] for s in list(raw_dc.detector_sources)]\n",
" modules = [int(s[:s.index(\"C\")]) for s in strings]\n",
" modules.sort()\n",
" karabo_da = [\"AGIPD{:02d}\".format(i) for i in modules]\n",
"else:\n",
" modules = [int(x[-2:]) for x in karabo_da]\n",
"\n",
"\n",
"print(\"Modules: {}\".format(modules))\n",
"print(f\"Detector in use is {karabo_id}\")"
]
},
{
"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": [
"# Evaluate creation time\n",
"creation_time = calcat_creation_time(in_folder, first_run, creation_time)\n",
"offset = parser.parse(creation_date_offset)\n",
"delta = timedelta(hours=offset.hour, minutes=offset.minute, seconds=offset.second)\n",
"creation_time += delta\n",
"print(f\"Creation time: {creation_time}\")\n",
"\n",
"# Read AGIPD parameter conditions.\n",
"if acq_rate == -1.:\n",
" acq_rate = agipd_cond.get_acq_rate()\n",
"if mem_cells == -1:\n",
" mem_cells = agipd_cond.get_num_cells()\n",
"# TODO: look for alternative for passing creation_time\n",
"if gain_setting == -1:\n",
" gain_setting = agipd_cond.get_gain_setting(creation_time)\n",
"if bias_voltage == -1:\n",
" bias_voltage = agipd_cond.get_bias_voltage(karabo_id_control)\n",
"if integration_time == -1:\n",
" integration_time = agipd_cond.get_integration_time()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(\"Operating conditions are:\")\n",
"print(f\"Bias voltage: {bias_voltage} V\")\n",
"print(f\"Memory cells: {mem_cells}\")\n",
"print(f\"Acquisition rate: {acq_rate} MHz\")\n",
"print(f\"Gain setting: {gain_setting}\")\n",
"print(f\"Integration time: {integration_time}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def count_min_bursts(in_folder, channel, run):\n",
" \"\"\"Calculate maximum number of trains in the provided runs.\"\"\"\n",
" \n",
" bursts = []\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",
" 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",
" return bursts\n",
"\n",
"trains = []\n",
"\n",
"for module in modules:\n",
" partial_check = partial(count_min_bursts, in_folder+\"/r{:04d}/\", module)\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 in runs: ',trains)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def check_ASIC(in_folder, runs, trains, channel, tresh):\n",
" \"\"\"Check which set of runs to use for a given ASIC\"\"\"\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*mem_cells, 0, :, :]\n",
" d = d.values.reshape((d.shape[0] // mem_cells, mem_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 is used to check increase of signal, \n",
" #if it is low, data from other file is checked\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",
" if (rise_val < tresh) & (rise_val > -30.):\n",
" ASIC_code.append(2)\n",
" print(f\"Signal rise low ({rise_val})\")\n",
" else:\n",
" ASIC_code.append(1)\n",
" print(f\"Signal rise ok ({rise_val})\")\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",
" if (rise_val < tresh) & (rise_val > -30.):\n",
" ASIC_code.append(2)\n",
" print(f\"Signal rise low ({rise_val})\")\n",
" else:\n",
" ASIC_code.append(1)\n",
" print(f\"Signal rise ok ({rise_val})\")\n",
"\n",
" return ASIC_code\n",
"\n",
"\n",
"\n",
"def defineInitParams(m_num):\n",
" \"\"\"Specify the ASIC indexing and runs to be used within data array in OrderedDict.\"\"\"\n",
" \n",
" channel = modules[m_num]\n",
" runs_dict = OrderedDict()\n",
" \n",
" instrument_source = karabo_id+'/DET/{}CH0:xtdf'.format(channel)\n",
" print(instrument_source)\n",
" \n",
" run_lbl = [0, runs1, runs2, runs3, runs4]\n",
" \n",
" for idx, a_c in enumerate(ASIC_code_list[m_num]):\n",
" \n",
" print(f'ASIC {idx} - 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": [
"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",
" substitute = np.where(ASIC_code2 == 2)[0]\n",
" ASIC_code[substitute] = 3\n",
" if 3 in ASIC_code:\n",
" print('\\n')\n",
" ASIC_code3 = check_ASIC(in_folder+\"/r{:04d}/\", runs3, trains, module, 100.)\n",
" ASIC_code3 = np.array(ASIC_code3)\n",
" ASIC_code3[np.where(ASIC_code == 1)[0]] = 1\n",
" ASIC_code3[np.where(ASIC_code == 2)[0]] = 5\n",
" if 2 in ASIC_code3:\n",
" substitute = np.where(ASIC_code3 == 2)[0]\n",
" ASIC_code[substitute] = 4\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": [
"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",
" \"\"\"Merge runs to create complete image\"\"\"\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",
" } \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(\"Processing Module: {}, ASIC: {}, Index: {}\".format(module, asic, index))\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",
" 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",
" \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": [
"start = datetime.now()\n",
"\n",
"run_lbl = [0, runs1, runs2, runs3, runs4]\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",
"\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",
" 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",
" f.flush()\n",
" \n",
" end = datetime.now() - start\n",
" print(f'Duration of merging for module {channel}: ', end)"
]
},
{
"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 = 4\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 single 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 = 4\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 = 4\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()"
]
}
],
"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
}