{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ePIX Data Correction ##\n",
"\n",
"Authors: Q. Tian S. Hauf, Version 1.0\n",
"\n",
"The following notebook provides Offset correction of images acquired with the ePix100 detector."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cluster_profile = \"noDB\" # ipcluster profile to use\n",
"in_folder = \"/gpfs/exfel/exp/CALLAB/202031/p900113/raw\" # input folder, required\n",
"out_folder = \"\" # output folder, required\n",
"sequences = [-1] # sequences to correct, set to -1 for all, range allowed\n",
"run = 9988 # which run to read data from, required\n",
"\n",
"karabo_id = \"MID_EXP_EPIX-1\" # karabo karabo_id\n",
"karabo_da = \"EPIX01\" # data aggregators\n",
"db_module = \"ePix100_M15\" # module id in the database\n",
"receiver_id = \"RECEIVER\" # inset for receiver devices\n",
"path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data\n",
"h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data/image' # path in the HDF5 file to images\n",
"h5path_t = 'data.backTemp' # path to find temperature at # TODO: modify in calibration_configurations.\n",
"h5path_cntrl = '/CONTROL/{}/DET' # path to control data # TODO: remove from calibration_configurations.\n",
"\n",
"use_dir_creation_date = True # date constants injected before directory creation time\n",
"cal_db_interface = \"tcp://max-exfl016:8015#8025\" # calibration DB interface to use\n",
"cal_db_timeout = 300000 # timeout on caldb requests\n",
"\n",
"cpuCores = 4 # Specifies the number of running cpu cores\n",
"chunk_size_idim = 1 # H5 chunking size of output data\n",
"overwrite = True # overwrite output folder\n",
"limit_images = 0 # process only first N images, 0 - process all\n",
"sequences_per_node = 1 # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel\n",
"\n",
"bias_voltage = 200 # bias voltage\n",
"in_vacuum = False # detector operated in vacuum\n",
"fix_temperature = 290. # fix temperature to this value\n",
"gain_photon_energy = 9.0 # Photon energy used for gain calibration\n",
"photon_energy = 8.0 # Photon energy to calibrate in number of photons, 0 for calibration in keV\n",
"\n",
"relative_gain = False # Apply relative gain correction.\n",
"\n",
"split_evt_primary_threshold = 7. # primary threshold for split event correction\n",
"split_evt_secondary_threshold = 5. # secondary threshold for split event correction\n",
"split_evt_mip_threshold = 1000. # minimum ionizing particle threshold\n",
"\n",
"\n",
"def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da):\n",
" from xfel_calibrate.calibrate import balance_sequences as bs\n",
" return bs(in_folder, run, sequences, sequences_per_node, karabo_da)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import tabulate\n",
"import warnings\n",
"\n",
"import h5py\n",
"import numpy as np\n",
"from IPython.display import Latex, display\n",
"from extra_data import RunDirectory\n",
"from pathlib import Path\n",
"\n",
"import XFELDetAna.xfelprofiler as xprof\n",
"from XFELDetAna import xfelpyanatools as xana\n",
"from XFELDetAna import xfelpycaltools as xcal\n",
"from XFELDetAna.plotting.util import prettyPlotting\n",
"from XFELDetAna.util import env\n",
"from cal_tools.tools import (\n",
" get_constant_from_db,\n",
" get_dir_creation_date,\n",
")\n",
"from iCalibrationDB import (\n",
" Conditions,\n",
" Constants,\n",
")\n",
"\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"prettyPlotting = True\n",
"\n",
"profiler = xprof.Profiler()\n",
"profiler.disable()\n",
"env.iprofile = cluster_profile\n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# TODO: expose to first cell after fixing common mode correction.\n",
"common_mode = False # Apply common mode correction.\n",
"\n",
"h5path = h5path.format(karabo_id, receiver_id)\n",
"plot_unit = 'ADU'\n",
"\n",
"if relative_gain:\n",
" plot_unit = 'keV'\n",
" if photon_energy > 0:\n",
" plot_unit = '$\\gamma$'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x = 708 # rows of the ePix100\n",
"y = 768 # columns of the ePix100\n",
"\n",
"in_folder = Path(in_folder)\n",
"ped_dir = in_folder / f\"r{run:04d}\"\n",
"run_dir = RunDirectory(ped_dir)\n",
"fp_name = path_template.format(run, karabo_da)\n",
"\n",
"print(f\"Reading from: {ped_dir / fp_name}\")\n",
"print(f\"Run is: {run}\")\n",
"print(f\"HDF5 path: {h5path}\")\n",
"print(f\"Data is output to: {out_folder}\")\n",
"\n",
"creation_time = None\n",
"if use_dir_creation_date:\n",
" creation_time = get_dir_creation_date(in_folder, run)\n",
"if creation_time:\n",
" print(f\"Using {creation_time.isoformat()} as creation time\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sensorSize = [x, y]\n",
"chunkSize = 100 # Number of images to read per chunk\n",
"# Sensor area will be analysed according to blocksize\n",
"blockSize = [sensorSize[0]//2, sensorSize[1]//2]\n",
"xcal.defaultBlockSize = blockSize\n",
"memoryCells = 1 # ePIX has no memory cells\n",
"run_parallel = True\n",
"\n",
"# Read slow data from the first available sequence file.\n",
"filename = ped_dir / fp_name.format(sequences[0] if sequences[0] != -1 else 0)\n",
"\n",
"integration_time = int(run_dir.get_array(\n",
" f\"{karabo_id}/DET/CONTROL\",\n",
" \"expTime.value\")[0])\n",
"temperature = np.mean(run_dir.get_array(\n",
" f\"{karabo_id}/DET/{receiver_id}:daqOutput\",\n",
" f\"{h5path_t}\").values) / 100.\n",
"\n",
"if fix_temperature != 0:\n",
" temperature_k = fix_temperature\n",
" print(\"Temperature is fixed!\")\n",
"else:\n",
" temperature_k = temperature + 273.15\n",
"\n",
"print(f\"Bias voltage is {bias_voltage} V\")\n",
"print(f\"Detector integration time is set to {integration_time}\")\n",
"print(\n",
" f\"Mean temperature was {temperature:0.2f} °C \"\n",
" f\"/ {temperature_k:0.2f} K at beginning of the run.\"\n",
")\n",
"print(f\"Operated in vacuum: {in_vacuum} \")\n",
"\n",
"out_folder = Path(out_folder)\n",
"if out_folder.is_dir() and not overwrite:\n",
" raise AttributeError(\"Output path exists! Exiting\")\n",
"out_folder.mkdir(parents=True, exist_ok=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Glob the right *.h5 fast data files.\n",
"seq_files = sorted(ped_dir.glob(f\"*{karabo_da}*.h5\"))\n",
"\n",
"# If a set of sequences requested to correct,\n",
"# adapt seq_files list.\n",
"if sequences != [-1]:\n",
" seq_files = [f for f in seq_files\n",
" if any(f.match(f\"*-S{s:05d}.h5\") for s in sequences)]\n",
"\n",
"print(f\"Processing a total of {len(seq_files)} sequence files\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Table of sequence files to process\n",
"table = [(k, f) for k, f in enumerate(seq_files)]\n",
"\n",
"if len(table):\n",
" md = display(Latex(tabulate.tabulate(\n",
" table,\n",
" tablefmt='latex',\n",
" headers=[\"#\", \"file\"]\n",
" )))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As a first step, dark maps have to be loaded."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"temp_limits = 5.\n",
"\n",
"cond_dict = {\n",
" \"bias_voltage\": bias_voltage,\n",
" \"integration_time\": integration_time,\n",
" \"temperature\": temperature_k,\n",
" \"in_vacuum\": in_vacuum,\n",
"}\n",
"\n",
"dark_condition = Conditions.Dark.ePix100(**cond_dict)\n",
"\n",
"# update conditions with illuminated conditins.\n",
"cond_dict.update({\n",
" \"photon_energy\": gain_photon_energy\n",
" })\n",
"\n",
"illum_condition = Conditions.Illuminated.ePix100(**cond_dict)\n",
"\n",
"const_cond = {\n",
" \"Offset\": dark_condition,\n",
" \"Noise\": dark_condition,\n",
" \"RelativeGain\": illum_condition,\n",
"}\n",
"\n",
"const_data = dict()\n",
"\n",
"for cname, condition in const_cond.items():\n",
" if cname == \"RelativeGain\" and not relative_gain:\n",
" continue\n",
" # TODO: Fix this logic.\n",
" for parm in condition.parameters:\n",
" if parm.name == \"Sensor Temperature\":\n",
" parm.lower_deviation = temp_limits\n",
" parm.upper_deviation = temp_limits\n",
"\n",
" const_data[cname] = get_constant_from_db(\n",
" karabo_id=karabo_id,\n",
" karabo_da=karabo_da,\n",
" constant=getattr(Constants.ePix100, cname)(),\n",
" condition=condition,\n",
" empty_constant=None,\n",
" cal_db_interface=cal_db_interface,\n",
" creation_time=creation_time,\n",
" print_once=2,\n",
" timeout=cal_db_timeout\n",
" )\n",
"\n",
"if relative_gain and const_data[\"RelativeGain\"] is None:\n",
" print(\n",
" \"WARNING: RelativeGain map is requested, but not found./n\"\n",
" \"No gain correction will be applied\"\n",
" )\n",
" relative_gain = False\n",
" plot_unit = 'ADU'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ************************Calculators************************ #\n",
"offsetCorrection = xcal.OffsetCorrection(\n",
" sensorSize,\n",
" const_data[\"Offset\"],\n",
" nCells=memoryCells,\n",
" cores=cpuCores,\n",
" gains=None,\n",
" blockSize=blockSize,\n",
" parallel=run_parallel\n",
")\n",
"\n",
"if relative_gain:\n",
" gainCorrection = xcal.RelativeGainCorrection(\n",
" sensorSize,\n",
" 1./const_data[\"RelativeGain\"][..., None],\n",
" nCells=memoryCells,\n",
" parallel=run_parallel,\n",
" cores=cpuCores,\n",
" blockSize=blockSize,\n",
" gains=None,\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# *****************Histogram Calculators****************** #\n",
"histCalOffsetCor = xcal.HistogramCalculator(\n",
" sensorSize,\n",
" bins=1050,\n",
" range=[-50, 1000],\n",
" parallel=run_parallel,\n",
" nCells=memoryCells,\n",
" cores=cpuCores,\n",
" blockSize=blockSize\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Applying corrections"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"histCalOffsetCor.debug()\n",
"offsetCorrection.debug()\n",
"if relative_gain:\n",
" gainCorrection.debug()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ************************Calculators************************ #\n",
"if common_mode:\n",
" commonModeBlockSize = [x//2, y//2]\n",
" commonModeAxisR = 'row'\n",
" cmCorrection = xcal.CommonModeCorrection(\n",
" [x, y],\n",
" commonModeBlockSize,\n",
" commonModeAxisR,\n",
" nCells=memoryCells,\n",
" noiseMap=const_data[\"Noise\"],\n",
" runParallel=run_parallel,\n",
" stats=True,\n",
" )\n",
"\n",
" histCalCMCor = xcal.HistogramCalculator(\n",
" sensorSize,\n",
" bins=1050,\n",
" range=[-50, 1000],\n",
" parallel=run_parallel,\n",
" nCells=memoryCells,\n",
" cores=cpuCores,\n",
" blockSize=blockSize,\n",
" )\n",
"\n",
"patternClassifier = xcal.PatternClassifier(\n",
" [x, y],\n",
" const_data[\"Noise\"],\n",
" split_evt_primary_threshold,\n",
" split_evt_secondary_threshold,\n",
" split_evt_mip_threshold,\n",
" tagFirstSingles=0,\n",
" nCells=memoryCells,\n",
" cores=cpuCores,\n",
" allowElongated=False,\n",
" blockSize=[x, y],\n",
" runParallel=run_parallel,\n",
")\n",
"\n",
"histCalSECor = xcal.HistogramCalculator(\n",
" sensorSize,\n",
" bins=1050,\n",
" range=[-50, 1000],\n",
" parallel=run_parallel,\n",
" nCells=memoryCells,\n",
" cores=cpuCores,\n",
" blockSize=blockSize,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"if common_mode:\n",
" cmCorrection.debug()\n",
" histCalCMCor.debug()\n",
"patternClassifier.debug()\n",
"histCalSECor.debug()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def copy_and_sanitize_non_cal_data(\n",
" infile: h5py,\n",
" outfile: h5py,\n",
" h5base: str\n",
"):\n",
" with h5py.File(f, \"r\") as infile:\n",
" \"\"\" Copy and sanitize data in `infile`\n",
" that is not touched by `correctEPIX`. \"\"\"\n",
"\n",
" if h5base.startswith(\"/\"):\n",
" h5base = h5base[1:]\n",
" dont_copy = [h5base+\"/pixels\"]\n",
"\n",
" def visitor(k, item):\n",
" if k not in dont_copy:\n",
" if isinstance(item, h5py.Group):\n",
" outfile.create_group(k)\n",
" elif isinstance(item, h5py.Dataset):\n",
" group = str(k).split(\"/\")\n",
" group = \"/\".join(group[:-1])\n",
" infile.copy(k, outfile[group])\n",
"\n",
" infile.visititems(visitor)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for f in seq_files:\n",
" data = None\n",
" out_file = out_folder / f.name.replace(\"RAW\", \"CORR\")\n",
" with h5py.File(out_file, \"w\") as ofile:\n",
" try:\n",
" copy_and_sanitize_non_cal_data(f, ofile, h5path)\n",
" data = run_dir.get_array(\n",
" f\"{karabo_id}/DET/{receiver_id}:daqOutput\",\n",
" \"data.image.pixels\").values\n",
" data = np.compress(np.any(data > 0, axis=(1, 2)), data, axis=0)\n",
" if limit_images > 0:\n",
" data = data[:limit_images, ...]\n",
"\n",
" oshape = data.shape\n",
" data = np.moveaxis(data, 0, 2)\n",
" ddset = ofile.create_dataset(\n",
" h5path+\"/pixels\",\n",
" oshape,\n",
" chunks=(chunk_size_idim, oshape[1], oshape[2]),\n",
" dtype=np.float32)\n",
"\n",
" # Offset correction.\n",
" data = offsetCorrection.correct(data.astype(np.float32))\n",
"\n",
" # relative gain correction.\n",
" if relative_gain:\n",
" data = gainCorrection.correct(data.astype(np.float32))\n",
" if photon_energy > 0:\n",
" data /= photon_energy\n",
"\n",
" histCalOffsetCor.fill(data)\n",
" ddset[...] = np.moveaxis(data, 2, 0)\n",
"\n",
" # Common mode correction\n",
" # TODO: Fix conflict between common mode and relative gain.\n",
"\n",
" \"\"\"The gain correction is currently applying an absolute correction\n",
" (not a relative correction as the implied by the name);\n",
" it changes the scale (the unit of measurement) of the data from ADU\n",
" to either keV or n_of_photons. But the common mode correction\n",
" relies on comparing data with the noise map, which is still in ADU.\n",
"\n",
" The best solution is probably to do a relative gain correction first\n",
" (correct) and apply the global absolute gain to the data at the end,\n",
" after common mode and clustering.\n",
" \"\"\"\n",
" if common_mode:\n",
" ddsetcm = ofile.create_dataset(\n",
" h5path+\"/pixels_cm\",\n",
" oshape,\n",
" chunks=(chunk_size_idim, oshape[1], oshape[2]),\n",
" dtype=np.float32)\n",
"\n",
" ddsetc = ofile.create_dataset(\n",
" h5path+\"/pixels_classified\",\n",
" oshape,\n",
" chunks=(chunk_size_idim, oshape[1], oshape[2]),\n",
" dtype=np.float32, compression=\"gzip\")\n",
"\n",
" ddsetp = ofile.create_dataset(\n",
" h5path+\"/patterns\",\n",
" oshape,\n",
" chunks=(chunk_size_idim, oshape[1], oshape[2]),\n",
" dtype=np.int32, compression=\"gzip\")\n",
"\n",
" # row common mode correction.\n",
" data = cmCorrection.correct(data)\n",
" histCalCMCor.fill(data)\n",
" ddsetcm[...] = np.moveaxis(data, 2, 0)\n",
"\n",
" data, patterns = patternClassifier.classify(data)\n",
"\n",
" data[data < (split_evt_primary_threshold*const_data[\"Noise\"])] = 0 # noqa\n",
" ddsetc[...] = np.moveaxis(data, 2, 0)\n",
" ddsetp[...] = np.moveaxis(patterns, 2, 0)\n",
"\n",
" data[patterns != 100] = np.nan\n",
" histCalSECor.fill(data)\n",
" except Exception as e:\n",
" print(f\"ERROR applying common mode correction for {f}: {e}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ho, eo, co, so = histCalOffsetCor.get()\n",
"\n",
"d = [{\n",
" 'x': co,\n",
" 'y': ho,\n",
" 'y_err': np.sqrt(ho[:]),\n",
" 'drawstyle': 'steps-mid',\n",
" 'errorstyle': 'bars',\n",
" 'errorcoarsing': 2,\n",
" 'label': 'Offset corr.'\n",
"}]\n",
"\n",
"if common_mode:\n",
" ho, eo, co, so = histCalCMCor.get()\n",
" d.append({\n",
" 'x': co,\n",
" 'y': ho,\n",
" 'y_err': np.sqrt(ho[:]),\n",
" 'drawstyle': 'steps-mid',\n",
" 'errorstyle': 'bars',\n",
" 'errorcoarsing': 2,\n",
" 'label': 'CM corr.'\n",
" })\n",
"\n",
" ho, eo, co, so = histCalSECor.get()\n",
" d.append({\n",
" 'x': co,\n",
" 'y': ho,\n",
" 'y_err': np.sqrt(ho[:]),\n",
" 'drawstyle': 'steps-mid',\n",
" 'errorstyle': 'bars',\n",
" 'errorcoarsing': 2,\n",
" 'label': 'Single split events'\n",
" })\n",
"\n",
"\n",
"fig = xana.simplePlot(\n",
" d, aspect=1, x_label=f'Energy({plot_unit})',\n",
" y_label='Number of occurrences', figsize='2col',\n",
" y_log=True, x_range=(-50, 500),\n",
" legend='top-center-frame-2col'\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mean Image of last Sequence ##"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = xana.heatmapPlot(\n",
" np.nanmedian(data, axis=2),\n",
" x_label='Columns', y_label='Rows',\n",
" lut_label=f'Signal ({plot_unit})',\n",
" x_range=(0, y),\n",
" y_range=(0, x),\n",
" vmin=-50, vmax=50)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Single Shot of last Sequence ##"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = xana.heatmapPlot(\n",
" data[..., 0],\n",
" x_label='Columns', y_label='Rows',\n",
" lut_label=f'Signal ({plot_unit})',\n",
" x_range=(0, y),\n",
" y_range=(0, x),\n",
" vmin=-50, vmax=50)"
]
}
],
"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"
},
"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": 1
}