diff --git a/notebooks/ePix100/Correction_ePix100_NBC.ipynb b/notebooks/ePix100/Correction_ePix100_NBC.ipynb
index bbfc6eae0920f0c4d937945c9d5c47d5f3713523..5df60c846b275c190a4352d5dcd1e8caec6ffa7c 100644
--- a/notebooks/ePix100/Correction_ePix100_NBC.ipynb
+++ b/notebooks/ePix100/Correction_ePix100_NBC.ipynb
@@ -17,7 +17,7 @@
"metadata": {},
"outputs": [],
"source": [
- "cluster_profile = \"noDB\" # ipcluster profile to use\n",
+ "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",
@@ -49,12 +49,12 @@
"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",
- "common_mode = False # Apply common mode 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",
+ "\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)"
@@ -66,37 +66,27 @@
"metadata": {},
"outputs": [],
"source": [
- "import copy\n",
- "import os\n",
"import tabulate\n",
- "import time\n",
"import warnings\n",
"\n",
"import h5py\n",
"import numpy as np\n",
- "from IPython.display import HTML, Latex, Markdown, display\n",
+ "from IPython.display import Latex, display\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.detectors.fastccd import readerh5 as fastccdreaderh5\n",
"from XFELDetAna.plotting.util import prettyPlotting\n",
- "import XFELDetAna.xfelprofiler as xprof\n",
- "from XFELDetAna.xfelreaders import ChunkReader\n",
"from XFELDetAna.util import env\n",
- "\n",
"from cal_tools.tools import (\n",
" get_constant_from_db,\n",
" get_dir_creation_date,\n",
")\n",
"from iCalibrationDB import (\n",
" Conditions,\n",
- " ConstantMetaData,\n",
" Constants,\n",
- " Detectors,\n",
- " Versions,\n",
")\n",
- "from iCalibrationDB.detectors import DetectorTypes\n",
"\n",
"warnings.filterwarnings('ignore')\n",
"\n",
@@ -115,11 +105,14 @@
"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",
"h5path_t = h5path_t.format(karabo_id, receiver_id)\n",
"h5path_cntrl = h5path_cntrl.format(karabo_id)\n",
- "\n",
"plot_unit = 'ADU'\n",
+ "\n",
"if relative_gain:\n",
" plot_unit = 'keV'\n",
" if photon_energy > 0:\n",
@@ -144,13 +137,11 @@
"print(f\"HDF5 path: {h5path}\")\n",
"print(f\"Data is output to: {out_folder}\")\n",
"\n",
- "import datetime\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\") "
+ " print(f\"Using {creation_time.isoformat()} as creation time\")"
]
},
{
@@ -161,23 +152,27 @@
"source": [
"sensorSize = [x, y]\n",
"chunkSize = 100 # Number of images to read per chunk\n",
- "blockSize = [sensorSize[0]//2, sensorSize[1]//2] # Sensor area will be analysed according to blocksize\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",
+ "filename = ped_dir / fp_name.format(sequences[0] if sequences[0] != -1 else 0)\n",
"with h5py.File(filename, 'r') as f:\n",
" integration_time = int(f[f\"{h5path_cntrl}/CONTROL/expTime/value\"][0])\n",
- " temperature = np.mean(f[h5path_t])/100.\n",
+ " temperature = np.mean(f[h5path_t]) / 100.\n",
" temperature_k = temperature + 273.15\n",
" if fix_temperature != 0:\n",
" temperature_k = fix_temperature\n",
" print(\"Temperature is fixed!\")\n",
" print(f\"Bias voltage is {bias_voltage} V\")\n",
" print(f\"Detector integration time is set to {integration_time}\")\n",
- " print(f\"Mean temperature was {temperature:0.2f} °C / {temperature_k:0.2f} K at beginning of run\")\n",
+ " print(\n",
+ " f\"Mean temperature was {temperature:0.2f} °C \"\n",
+ " f\"/ {temperature_k:0.2f} K at beginning of run\"\n",
+ " )\n",
" print(f\"Operated in vacuum: {in_vacuum} \")\n",
"\n",
"out_folder = Path(out_folder)\n",
@@ -198,7 +193,8 @@
"# 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 if any(f.match(f\"*-S{s:05d}.h5\") for s in sequences)]\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\")"
]
@@ -212,8 +208,12 @@
"# 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(table, tablefmt='latex', headers=[\"#\", \"file\"]))) "
+ "if len(table):\n",
+ " md = display(Latex(tabulate.tabulate(\n",
+ " table,\n",
+ " tablefmt='latex',\n",
+ " headers=[\"#\", \"file\"]\n",
+ " )))"
]
},
{
@@ -238,7 +238,7 @@
" \"in_vacuum\": in_vacuum,\n",
"}\n",
"\n",
- "dark_condition = Conditions.Dark.ePix100(**cond_dict)\n",
+ "dark_condition = Conditions.Dark.ePix100(**cond_dict)\n",
"\n",
"# update conditions with illuminated conditins.\n",
"cond_dict.update({\n",
@@ -291,11 +291,11 @@
"metadata": {},
"outputs": [],
"source": [
- "#************************Calculators************************#\n",
+ "# ************************Calculators************************ #\n",
"offsetCorrection = xcal.OffsetCorrection(\n",
" sensorSize,\n",
" const_data[\"Offset\"],\n",
- " nCells = memoryCells,\n",
+ " nCells=memoryCells,\n",
" cores=cpuCores,\n",
" gains=None,\n",
" blockSize=blockSize,\n",
@@ -304,7 +304,7 @@
"\n",
"if relative_gain:\n",
" gainCorrection = xcal.RelativeGainCorrection(\n",
- " sensorSize, \n",
+ " sensorSize,\n",
" 1./const_data[\"RelativeGain\"][..., None],\n",
" nCells=memoryCells,\n",
" parallel=run_parallel,\n",
@@ -320,7 +320,7 @@
"metadata": {},
"outputs": [],
"source": [
- "#*****************Histogram Calculators******************#\n",
+ "# *****************Histogram Calculators****************** #\n",
"histCalOffsetCor = xcal.HistogramCalculator(\n",
" sensorSize,\n",
" bins=1050,\n",
@@ -357,7 +357,7 @@
"metadata": {},
"outputs": [],
"source": [
- "#************************Calculators************************#\n",
+ "# ************************Calculators************************ #\n",
"if common_mode:\n",
" commonModeBlockSize = [x//2, y//2]\n",
" commonModeAxisR = 'row'\n",
@@ -372,35 +372,35 @@
" )\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",
+ " sensorSize,\n",
+ " bins=1050,\n",
+ " range=[-50, 1000],\n",
+ " parallel=run_parallel,\n",
+ " nCells=memoryCells,\n",
+ " cores=cpuCores,\n",
+ " blockSize=blockSize,\n",
" )\n",
- " \n",
+ "\n",
"patternClassifier = xcal.PatternClassifier(\n",
- " [x, y], \n",
- " const_data[\"Noise\"], \n",
- " split_evt_primary_threshold, \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",
+ " 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",
+ " sensorSize,\n",
+ " bins=1050,\n",
" range=[-50, 1000],\n",
" parallel=run_parallel,\n",
- " nCells=memoryCells, \n",
+ " nCells=memoryCells,\n",
" cores=cpuCores,\n",
" blockSize=blockSize,\n",
")"
@@ -430,9 +430,9 @@
" outfile: h5py,\n",
" h5base: str\n",
"):\n",
- " \"\"\" Copy and sanitize data in `infile` that is not touched by `correctEPIX`\n",
- " \"\"\"\n",
- " \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",
@@ -445,7 +445,7 @@
" group = str(k).split(\"/\")\n",
" group = \"/\".join(group[:-1])\n",
" infile.copy(k, outfile[group])\n",
- " \n",
+ "\n",
" infile.visititems(visitor)"
]
},
@@ -462,9 +462,9 @@
" try:\n",
" copy_and_sanitize_non_cal_data(infile, ofile, h5path)\n",
" data = infile[h5path+\"/pixels\"][()]\n",
- " data = np.compress(np.any(data>0, axis=(1,2)), data, axis=0)\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",
+ " data = data[:limit_images, ...]\n",
"\n",
" oshape = data.shape\n",
" data = np.moveaxis(data, 0, 2)\n",
@@ -487,6 +487,18 @@
" 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",
@@ -507,13 +519,13 @@
" dtype=np.int32, compression=\"gzip\")\n",
"\n",
" # row common mode correction.\n",
- " data = cmCorrection.correct(data) \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\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",
@@ -529,7 +541,7 @@
"metadata": {},
"outputs": [],
"source": [
- "ho,eo,co,so = histCalOffsetCor.get()\n",
+ "ho, eo, co, so = histCalOffsetCor.get()\n",
"\n",
"d = [{\n",
" 'x': co,\n",
@@ -542,7 +554,7 @@
"}]\n",
"\n",
"if common_mode:\n",
- " ho,eo,co,so = histCalCMCor.get()\n",
+ " ho, eo, co, so = histCalCMCor.get()\n",
" d.append({\n",
" 'x': co,\n",
" 'y': ho,\n",
@@ -553,7 +565,7 @@
" 'label': 'CM corr.'\n",
" })\n",
"\n",
- " ho,eo,co,so = histCalSECor.get()\n",
+ " ho, eo, co, so = histCalSECor.get()\n",
" d.append({\n",
" 'x': co,\n",
" 'y': ho,\n",
@@ -568,7 +580,7 @@
"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",
+ " y_log=True, x_range=(-50, 500),\n",
" legend='top-center-frame-2col'\n",
")"
]
@@ -590,8 +602,8 @@
" 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",
+ " x_range=(0, y),\n",
+ " y_range=(0, x),\n",
" vmin=-50, vmax=50)"
]
},
@@ -609,11 +621,11 @@
"outputs": [],
"source": [
"fig = xana.heatmapPlot(\n",
- " data[...,0],\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",
+ " x_range=(0, y),\n",
+ " y_range=(0, x),\n",
" vmin=-50, vmax=50)"
]
}