diff --git a/notebooks/ePix100/Correction_ePix100_NBC.ipynb b/notebooks/ePix100/Correction_ePix100_NBC.ipynb index 5158373979d84e14ac643d0c5f61dbd6ccaab81c..2f6e0d90b5dc28443fc38b9bf14a361904810ef1 100644 --- a/notebooks/ePix100/Correction_ePix100_NBC.ipynb +++ b/notebooks/ePix100/Correction_ePix100_NBC.ipynb @@ -24,15 +24,15 @@ "metadata": {}, "outputs": [], "source": [ - "in_folder = \"/gpfs/exfel/exp/HED/202202/p003121/raw\" # input folder, required\n", + "in_folder = \"/gpfs/exfel/exp/MID/202301/p003346/raw\" # input folder, required\n", "out_folder = \"\" # output folder, required\n", "metadata_folder = \"\" # Directory containing calibration_metadata.yml when run by xfel-calibrate\n", "sequences = [-1] # sequences to correct, set to -1 for all, range allowed\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", - "run = 156 # which run to read data from, required\n", + "run = 55 # which run to read data from, required\n", "\n", "# Parameters for accessing the raw data.\n", - "karabo_id = \"HED_IA1_EPX100-1\" # karabo karabo_id\n", + "karabo_id = \"MID_EXP_EPIX-1\" # karabo karabo_id\n", "karabo_da = \"EPIX01\" # data aggregators\n", "db_module = \"\" # module id in the database\n", "receiver_template = \"RECEIVER\" # detector receiver template for accessing raw data files\n", @@ -92,14 +92,16 @@ "import matplotlib.pyplot as plt\n", "from IPython.display import Latex, Markdown, display\n", "from extra_data import RunDirectory, H5File\n", + "from extra_geom import Epix100Geometry\n", + "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", "from pathlib import Path\n", "\n", "import cal_tools.restful_config as rest_cfg\n", - "from XFELDetAna import xfelpyanatools as xana\n", "from XFELDetAna import xfelpycaltools as xcal\n", "from cal_tools.calcat_interface import EPIX100_CalibrationData\n", "from cal_tools.epix100 import epix100lib\n", "from cal_tools.files import DataFile\n", + "from cal_tools.restful_config import restful_config\n", "from cal_tools.tools import (\n", " calcat_creation_time,\n", " write_constants_fragment,\n", @@ -110,6 +112,7 @@ "\n", "prettyPlotting = True\n", "\n", + "\n", "%matplotlib inline" ] }, @@ -134,6 +137,8 @@ "metadata": {}, "outputs": [], "source": [ + "prop_str = in_folder[in_folder.find('/p')+1:in_folder.find('/p')+8]\n", + "\n", "in_folder = Path(in_folder)\n", "out_folder = Path(out_folder)\n", "\n", @@ -250,13 +255,12 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Retrieving calibration constants\n", "\n", - "As a first step, calibration constants have to be loaded." + "As a first step, dark maps have to be loaded." ] }, { @@ -321,35 +325,8 @@ "# Initializing some parameters.\n", "hscale = 1\n", "stats = True\n", - "hrange = np.array([-50, 1000])\n", - "nbins = hrange[1] - hrange[0]\n", - "commonModeBlockSize = [x//2, y//2]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "histCalOffsetCor = xcal.HistogramCalculator(\n", - " sensorSize,\n", - " bins=nbins,\n", - " range=hrange,\n", - " parallel=run_parallel,\n", - " nCells=memoryCells,\n", - " blockSize=blockSize\n", - ")\n", - "\n", - "# *****************Histogram Calculators****************** #\n", - "histCalCor = xcal.HistogramCalculator(\n", - " sensorSize,\n", - " bins=1050,\n", - " range=[-50, 1000],\n", - " parallel=run_parallel,\n", - " nCells=memoryCells,\n", - " blockSize=blockSize\n", - ")" + "bins = np.arange(-50,1000)\n", + "hist = {'O': 0} # dictionary to store histograms" ] }, { @@ -359,14 +336,9 @@ "outputs": [], "source": [ "if common_mode:\n", - " histCalCMCor = xcal.HistogramCalculator(\n", - " sensorSize,\n", - " bins=nbins,\n", - " range=hrange,\n", - " parallel=run_parallel,\n", - " nCells=memoryCells,\n", - " blockSize=blockSize,\n", - " )\n", + " \n", + " commonModeBlockSize = [x//2, y//8]\n", + " \n", " cmCorrectionB = xcal.CommonModeCorrection(\n", " shape=sensorSize,\n", " blockSize=commonModeBlockSize, \n", @@ -402,7 +374,9 @@ " stats=stats,\n", " minFrac=cm_min_frac,\n", " noiseSigma=cm_noise_sigma,\n", - " )" + " )\n", + " \n", + " hist['CM'] = 0" ] }, { @@ -412,12 +386,11 @@ "outputs": [], "source": [ "if relative_gain:\n", - " gain_cnst = np.median(const_data[\"RelativeGainEPix100\"])\n", - " hscale = gain_cnst\n", - " plot_unit = 'keV'\n", - " if photon_energy > 0:\n", - " plot_unit = '$\\gamma$'\n", - " hscale /= photon_energy\n", + " \n", + " # gain constant is given by the mode of the gain map \n", + " # because all bad pixels are masked using this value\n", + " _vals,_counts = np.unique(const_data[\"RelativeGainEPix100\"], return_counts=True)\n", + " gain_cnst = _vals[np.argmax(_counts)] \n", " \n", " gainCorrection = xcal.RelativeGainCorrection(\n", " sensorSize,\n", @@ -427,25 +400,16 @@ " blockSize=blockSize,\n", " gains=None,\n", " )\n", - "\n", - " histCalRelGainCor = xcal.HistogramCalculator(\n", - " sensorSize,\n", - " bins=nbins,\n", - " range=hrange,\n", - " parallel=run_parallel,\n", - " nCells=memoryCells,\n", - " blockSize=blockSize\n", - " )\n", + " \n", + " hist['RG'] = 0\n", "\n", " if absolute_gain:\n", - " histCalAbsGainCor = xcal.HistogramCalculator(\n", - " sensorSize,\n", - " bins=nbins,\n", - " range=hrange*hscale,\n", - " parallel=run_parallel,\n", - " nCells=memoryCells,\n", - " blockSize=blockSize\n", - " )" + " hscale = gain_cnst\n", + " plot_unit = 'keV'\n", + " if photon_energy > 0:\n", + " plot_unit = '$\\gamma$'\n", + " hscale /= photon_energy\n", + " hist['AG'] = 0" ] }, { @@ -467,30 +431,8 @@ " blockSize=[x, y],\n", " parallel=run_parallel,\n", " )\n", - " histCalCSCor = xcal.HistogramCalculator(\n", - " sensorSize,\n", - " bins=nbins,\n", - " range=hrange,\n", - " parallel=run_parallel,\n", - " nCells=memoryCells,\n", - " blockSize=blockSize,\n", - " )\n", - " histCalGainCorClusters = xcal.HistogramCalculator(\n", - " sensorSize,\n", - " bins=nbins,\n", - " range=hrange*hscale,\n", - " parallel=run_parallel,\n", - " nCells=memoryCells,\n", - " blockSize=blockSize\n", - " )\n", - " histCalGainCorSingles = xcal.HistogramCalculator(\n", - " sensorSize,\n", - " bins=nbins,\n", - " range=hrange*hscale,\n", - " parallel=run_parallel,\n", - " nCells=memoryCells,\n", - " blockSize=blockSize\n", - " )" + " hist['CS'] = 0\n", + " hist['S'] = 0" ] }, { @@ -511,11 +453,11 @@ " d = d[..., np.newaxis].astype(np.float32)\n", " d = np.compress(\n", " np.any(d > 0, axis=(0, 1)), d, axis=2)\n", - " \n", + "\n", " # Offset correction.\n", " d -= const_data[\"OffsetEPix100\"]\n", + " hist['O'] += np.histogram(d,bins=bins)[0]\n", "\n", - " histCalOffsetCor.fill(d)\n", " # Common Mode correction.\n", " if common_mode:\n", " # Block CM\n", @@ -524,12 +466,15 @@ " d = cmCorrectionR.correct(d)\n", " # COL CM\n", " d = cmCorrectionC.correct(d)\n", - " histCalCMCor.fill(d)\n", "\n", - " # relative gain correction.\n", + " hist['CM'] += np.histogram(d,bins=bins)[0]\n", + "\n", + "\n", + " # Relative gain correction.\n", " if relative_gain:\n", " d = gainCorrection.correct(d)\n", - " histCalRelGainCor.fill(d)\n", + "\n", + " hist['RG'] += np.histogram(d,bins=bins)[0]\n", "\n", " \"\"\"The gain correction is currently applying\n", " an absolute correction (not a relative correction\n", @@ -552,34 +497,31 @@ " data_clu[index, ...] = np.squeeze(d_clu)\n", " data_patterns[index, ...] = np.squeeze(patterns)\n", "\n", - " histCalCSCor.fill(d_clu)\n", + " hist['CS'] += np.histogram(d_clu,bins=bins)[0]\n", "\n", - " # absolute gain correction\n", + " d_sing = d_clu[patterns==100] # pattern 100 corresponds to single photons events\n", + " if len(d_sing):\n", + " hist['S'] += np.histogram(d_sing,bins=bins)[0]\n", + "\n", + " # Absolute gain correction\n", " # changes data from ADU to keV (or n. of photons)\n", " if absolute_gain:\n", "\n", " d = d * gain_cnst\n", " if photon_energy > 0:\n", " d /= photon_energy\n", - " histCalAbsGainCor.fill(d)\n", + "\n", + " hist['AG'] += np.histogram(d,bins=bins)[0]\n", "\n", " if pattern_classification:\n", " # Modify pattern classification.\n", " d_clu = d_clu * gain_cnst\n", - " \n", " if photon_energy > 0:\n", " d_clu /= photon_energy\n", "\n", " data_clu[index, ...] = np.squeeze(d_clu)\n", "\n", - " histCalGainCorClusters.fill(d_clu)\n", - " \n", - " d_sing = d_clu[patterns==100] # pattern 100 corresponds to single photons events\n", - " if len(d_sing):\n", - " histCalGainCorSingles.fill(d_sing)\n", - "\n", - " data[index, ...] = np.squeeze(d)\n", - " histCalCor.fill(d)" + " data[index, ...] = np.squeeze(d)" ] }, { @@ -712,80 +654,52 @@ " exit(0)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot Histograms" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bins_ADU = bins[:-1]+np.diff(bins)[0]/2\n", + "bins_keV = bins_ADU*hscale" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "ho, eo, co, so = histCalCor.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': 'Total corr.'\n", - "}]\n", - "\n", - "ho, eo, co, so = histCalOffsetCor.get()\n", - "\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': 'Offset corr.'\n", - "})\n", + "# Histogram in ADU\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", - "if relative_gain :\n", - " ho, eo, co, so = histCalRelGainCor.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': 'Relative gain corr.'\n", - " })\n", + "plt.figure(figsize=(12,8))\n", + "plt.plot(bins_ADU,hist['O'], label='Offset corr')\n", "\n", + "if common_mode:\n", + " plt.plot(bins_ADU,hist['CM'], label='CM corr')\n", + "if relative_gain:\n", + " plt.plot(bins_ADU,hist['RG'], label='Relative Gain corr')\n", "if pattern_classification:\n", - " ho, eo, co, so = histCalCSCor.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': 'Charge sharing corr.'\n", - " })\n", - "\n", - "fig = xana.simplePlot(\n", - " d, aspect=1, x_label=f'Energy (ADU)',\n", - " y_label='Number of occurrences', figsize='2col',\n", - " y_log=True, x_range=(-50, 500),\n", - " legend='top-center-frame-2col',\n", - ")\n", - "plt.title(f'run {run} - {karabo_da}')\n", - "plt.grid()" + " plt.plot(bins_ADU[bins_ADU>10],hist['CS'][bins_ADU>10], label='Charge Sharing corr')\n", + " if np.any(hist['S']):\n", + " plt.plot(bins_ADU,hist['S'], label='Singles')\n", + "\n", + "xtick_step = 50\n", + "plt.xlim(bins[0], bins[-1]+1)\n", + "plt.xticks(np.arange(bins[0],bins[-1]+2,xtick_step))\n", + "plt.xlabel('ADU',fontsize=12)\n", + "\n", + "plt.yscale('log')\n", + "plt.title(f'{karabo_id} | {prop_str}, r{run}', fontsize=14, fontweight='bold')\n", + "plt.legend(fontsize=12)\n", + "plt.grid(ls=':')" ] }, { @@ -794,54 +708,30 @@ "metadata": {}, "outputs": [], "source": [ - "if absolute_gain :\n", - " d = []\n", - " ho, eo, co, so = histCalAbsGainCor.get()\n", - " if co is not None: # avoid adding None array, if calculator is empty.\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': 'Absolute gain corr.'\n", - " })\n", + "# Histogram in keV/number of photons\n", + "\n", + "if absolute_gain:\n", + " colors = plt.rcParams['axes.prop_cycle'].by_key()['color']\n", + " plt.figure(figsize=(12,8))\n", "\n", + " if relative_gain:\n", + " plt.plot(bins_keV,hist['RG'], label='Absolute Gain corr', c=colors[2])\n", + " \n", " if pattern_classification:\n", - " ho, eo, co, so = histCalGainCorClusters.get()\n", - " if co is not None: # avoid adding None array, if calculator is empty.\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': 'Charge sharing corr.'\n", - " })\n", - " \n", - " ho, eo, co, so = histCalGainCorSingles.get()\n", - " if co is not None: # avoid adding None array, if calculator is empty.\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': 'Isolated photons (singles)'\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, \n", - " x_range=np.array((-50, 500))*hscale,\n", - " legend='top-center-frame-2col',\n", - " )\n", - " plt.grid()\n", - " plt.title(f'run {run} - {karabo_da}')" + " plt.plot(bins_keV[bins_keV>.5],hist['CS'][bins_keV>.5], label='Charge Sharing corr', c=colors[3])\n", + " if np.any(hist['S']):\n", + " plt.plot(bins_keV[bins_keV>.5],hist['S'][bins_keV>.5], label='Singles', c=colors[4])\n", + " \n", + " if photon_energy==0: # if keV instead of #photons\n", + " xtick_step = 5\n", + " plt.xlim(left=-2)\n", + " plt.xticks(np.arange(0,plt.gca().get_xlim()[1],xtick_step))\n", + " plt.xlabel(plot_unit,fontsize=12)\n", + "\n", + " plt.yscale('log')\n", + " plt.title(f'{karabo_id} | {prop_str}, r{run}', fontsize=14, fontweight='bold')\n", + " plt.legend(fontsize=12)\n", + " plt.grid(ls=':')" ] }, { @@ -857,15 +747,48 @@ "metadata": {}, "outputs": [], "source": [ - "step_timer.start()\n", - "fig = xana.heatmapPlot(\n", - " np.nanmedian(data, axis=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)\n", - "step_timer.done_step(f'Plotting mean image of {data.shape[0]} trains.')" + "geom = Epix100Geometry.from_relative_positions(top=[386.5, 364.5, 0.], bottom=[386.5, -12.5, 0.])\n", + "\n", + "if pattern_classification:\n", + " plt.subplots(1,2,figsize=(18,18)) if pattern_classification else plt.subplots(1,1,figsize=(9,9))\n", + " ax = plt.subplot(1,2,1)\n", + " ax.set_title(f'Before CS correction',fontsize=12,fontweight='bold');\n", + "else:\n", + " plt.subplots(1,1,figsize=(9,9))\n", + " ax = plt.subplot(1,1,1)\n", + " ax.set_title(f'{karabo_id} | {prop_str}, r{run} | Average of {data.shape[0]} trains',fontsize=12,fontweight='bold');\n", + " \n", + "# Average image before charge sharing corrcetion\n", + "divider = make_axes_locatable(ax)\n", + "cax = divider.append_axes('bottom', size='5%', pad=0.5)\n", + "\n", + "image = data.mean(axis=0)\n", + "vmin = max(image.mean()-2*image.std(),0)\n", + "vmax = image.mean()+3*image.std()\n", + "geom.plot_data(image,\n", + " ax=ax,\n", + " colorbar={'cax': cax, 'label': plot_unit, 'orientation': 'horizontal'},\n", + " origin='upper',\n", + " vmin=vmin,\n", + " vmax=vmax)\n", + "\n", + "# Average image after charge sharing corrcetion\n", + "if pattern_classification:\n", + "\n", + " ax = plt.subplot(1,2,2)\n", + " divider = make_axes_locatable(ax)\n", + " cax = divider.append_axes('bottom', size='5%', pad=0.5)\n", + "\n", + " image = data_clu.mean(axis=0)\n", + " geom.plot_data(image,\n", + " ax=ax,\n", + " colorbar={'cax': cax, 'label': plot_unit, 'orientation': 'horizontal'},\n", + " origin='upper',\n", + " vmin=vmin,\n", + " vmax=vmax)\n", + " ax.set_title(f'After CS correction',fontsize=12,fontweight='bold');\n", + "\n", + " plt.suptitle(f'{karabo_id} | {prop_str}, r{run} | Average of {data.shape[0]} trains',fontsize=14,fontweight='bold',y=.72);" ] }, { @@ -881,23 +804,56 @@ "metadata": {}, "outputs": [], "source": [ - "step_timer.start()\n", - "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)\n", - "step_timer.done_step(f'Plotting single shot of corrected data.')" + "train_idx = -1\n", + "\n", + "if pattern_classification:\n", + " plt.subplots(1,2,figsize=(18,18)) if pattern_classification else plt.subplots(1,1,figsize=(9,9))\n", + " ax = plt.subplot(1,2,1)\n", + " ax.set_title(f'Before CS correction',fontsize=12,fontweight='bold');\n", + "else:\n", + " plt.subplots(1,1,figsize=(9,9))\n", + " ax = plt.subplot(1,1,1)\n", + " ax.set_title(f'{karabo_id} | {prop_str}, r{run} | Single frame',fontsize=12,fontweight='bold');\n", + " \n", + "# Average image before charge sharing corrcetion\n", + "divider = make_axes_locatable(ax)\n", + "cax = divider.append_axes('bottom', size='5%', pad=0.5)\n", + "\n", + "image = data[train_idx]\n", + "vmin = max(image.mean()-2*image.std(),0)\n", + "vmax = image.mean()+3*image.std()\n", + "geom.plot_data(image,\n", + " ax=ax,\n", + " colorbar={'cax': cax, 'label': plot_unit, 'orientation': 'horizontal'},\n", + " origin='upper',\n", + " vmin=vmin,\n", + " vmax=vmax)\n", + "\n", + "# Average image after charge sharing corrcetion\n", + "if pattern_classification:\n", + "\n", + " ax = plt.subplot(1,2,2)\n", + " divider = make_axes_locatable(ax)\n", + " cax = divider.append_axes('bottom', size='5%', pad=0.5)\n", + "\n", + " image = data_clu[train_idx]\n", + " geom.plot_data(image,\n", + " ax=ax,\n", + " colorbar={'cax': cax, 'label': plot_unit, 'orientation': 'horizontal'},\n", + " origin='upper',\n", + " vmin=vmin,\n", + " vmax=vmax)\n", + " ax.set_title(f'After CS correction',fontsize=12,fontweight='bold');\n", + "\n", + " plt.suptitle(f'{karabo_id} | {prop_str}, r{run} | Single frame',fontsize=14,fontweight='bold',y=.72);" ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "cal_venv", "language": "python", - "name": "python3" + "name": "cal_venv" }, "language_info": { "codemirror_mode": { @@ -909,7 +865,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.12" + "version": "3.8.11" }, "latex_envs": { "LaTeX_envs_menu_present": true,