diff --git a/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb b/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb index c4431088e4a5e693f9e89cd3464c8183471240b5..7ca70a974be7c05d321b7e5f994567414c28590b 100644 --- a/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb +++ b/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb @@ -105,6 +105,7 @@ "# Plotting parameters\n", "skip_plots = False # exit after writing corrected files and metadata\n", "cell_id_preview = 1 # cell Id used for preview in single-shot plots\n", + "cmap = \"viridis\" # matplolib.colormap for almost all heatmap. Other options ['plasma', 'inferno', 'magma', 'cividis', 'jet', ...]\n", "\n", "# Parallelization parameters\n", "chunk_size = 1000 # Size of chunk for image-wise correction\n", @@ -114,6 +115,7 @@ "max_nodes = 8 # Maximum number of SLURM jobs to split correction work into\n", "max_tasks_per_worker = 1 # the number of tasks a correction pool worker process can complete before it will exit and be replaced with a fresh worker process. Leave as -1 to keep worker alive as long as pool.\n", "\n", + "\n", "def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da, max_nodes):\n", " from xfel_calibrate.calibrate import balance_sequences as bs\n", " return bs(in_folder, run, sequences, sequences_per_node, karabo_da, max_nodes=max_nodes)" @@ -909,24 +911,7 @@ "metadata": {}, "outputs": [], "source": [ - "def do_3d_plot(data, edges, x_axis, y_axis):\n", - " fig = plt.figure(figsize=(10, 10))\n", - " ax = fig.gca(projection='3d')\n", - "\n", - " # Make data.\n", - " X = edges[0][:-1]\n", - " Y = edges[1][:-1]\n", - " X, Y = np.meshgrid(X, Y)\n", - " Z = data.T\n", - "\n", - " # Plot the surface.\n", - " ax.plot_surface(X, Y, Z, cmap=colormap.coolwarm, linewidth=0, antialiased=False)\n", - " ax.set_xlabel(x_axis)\n", - " ax.set_ylabel(y_axis)\n", - " ax.set_zlabel(\"Counts\")\n", - "\n", - "\n", - "def do_2d_plot(data, edges, y_axis, x_axis):\n", + "def do_2d_plot(data, edges, y_axis, x_axis, title=\"\"):\n", " fig = plt.figure(figsize=(10, 10))\n", " ax = fig.add_subplot(111)\n", " extent = [np.min(edges[1]), np.max(edges[1]),\n", @@ -935,6 +920,7 @@ " norm=LogNorm(vmin=1, vmax=max(10, np.max(data))))\n", " ax.set_xlabel(x_axis)\n", " ax.set_ylabel(y_axis)\n", + " ax.set_title(title)\n", " cb = fig.colorbar(im)\n", " cb.set_label(\"Counts\")" ] @@ -1031,12 +1017,19 @@ "metadata": {}, "outputs": [], "source": [ - "hist, bins_x, bins_y = calgs.histogram2d(raw[:,0,...].flatten().astype(np.float32),\n", - " raw[:,1,...].flatten().astype(np.float32),\n", - " bins=(100, 100),\n", - " range=[[4000, 8192], [4000, 8192]])\n", - "do_2d_plot(hist, (bins_x, bins_y), \"Signal (ADU)\", \"Analogue gain (ADU)\")\n", - "do_3d_plot(hist, (bins_x, bins_y), \"Signal (ADU)\", \"Analogue gain (ADU)\")" + "raw_float = raw.astype(np.float32)\n", + "signal = raw[:, 0, ...]\n", + "gain = raw[:, 1, ...]\n", + "hist, bins_x, bins_y = calgs.histogram2d(\n", + " signal.flatten().astype(np.float32),\n", + " gain.flatten().astype(np.float32),\n", + " bins=(100, 100),\n", + " range=[\n", + " np.percentile(signal, [0.02, 99.8]),\n", + " np.percentile(gain, [0.02, 99.8]),\n", + " ],\n", + " )\n", + "do_2d_plot(hist, (bins_x, bins_y), \"Signal (ADU)\", \"Analogue gain (ADU)\")" ] }, { @@ -1054,9 +1047,16 @@ "metadata": {}, "outputs": [], "source": [ - "hist, bins_x, bins_y = calgs.histogram2d(corrected.flatten().astype(np.float32),\n", - " gains.flatten().astype(np.float32), bins=(100, 3),\n", - " range=[[-50, 8192], [0, 3]])\n", + "vmin, vmax = np.nanmin(corrected), np.nanmax(corrected)\n", + "hist, bins_x, bins_y = calgs.histogram2d(\n", + " corrected.flatten().astype(np.float32),\n", + " gains.flatten().astype(np.float32), bins=(100, 3),\n", + " range=[\n", + " # The range boundaries and decided by DET expert.\n", + " [max(vmin, -50), min(vmax, 8192)],\n", + " [0, 3]\n", + " ],\n", + " )\n", "do_2d_plot(hist, (bins_x, bins_y), \"Signal (ADU)\", \"Gain bit value\")" ] }, @@ -1089,26 +1089,33 @@ "source": [ "pulse_range = [np.min(pulseId[pulseId>=0]), np.max(pulseId[pulseId>=0])]\n", "\n", + "\n", + "def clamp(value, min_value, max_value):\n", + " return max(min_value, min(value, max_value))\n", + "\n", + "\n", "# Modify pulse_range, if only one pulse is selected.\n", "if pulse_range[0] == pulse_range[1]:\n", " pulse_range = [0, pulse_range[1]+int(acq_rate)]\n", "\n", "mean_data = np.nanmean(corrected, axis=(2, 3))\n", - "hist, bins_x, bins_y = calgs.histogram2d(mean_data.flatten().astype(np.float32),\n", - " pulseId.flatten().astype(np.float32),\n", - " bins=(100, int(pulse_range[1])),\n", - " range=[[-50, 1000], pulse_range])\n", - "\n", - "do_2d_plot(hist, (bins_x, bins_y), \"Signal (ADU)\", \"Pulse id\")\n", - "do_3d_plot(hist, (bins_x, bins_y), \"Signal (ADU)\", \"Pulse id\")\n", - "\n", - "hist, bins_x, bins_y = calgs.histogram2d(mean_data.flatten().astype(np.float32),\n", - " pulseId.flatten().astype(np.float32),\n", - " bins=(100, int(pulse_range[1])),\n", - " range=[[-50, 200000], pulse_range])\n", + "vmin, vmax = mean_data.min(), mean_data.max()\n", + "hist, bins_x, bins_y = calgs.histogram2d(\n", + " mean_data.flatten().astype(np.float32),\n", + " pulseId.flatten().astype(np.float32),\n", + " bins=(100, int(pulse_range[1])),\n", + " range=[[clamp(vmin, -50, -0.2), min(vmax, 1000)], pulse_range],\n", + ")\n", + "do_2d_plot(hist, (bins_x, bins_y), \"Signal (ADU)\", \"Pulse id\", title=\"Signal-Pulse ID\")\n", "\n", - "do_2d_plot(hist, (bins_x, bins_y), \"Signal (ADU)\", \"Pulse id\")\n", - "do_3d_plot(hist, (bins_x, bins_y), \"Signal (ADU)\", \"Pulse id\")" + "if vmax > 1000: # a zoom out plot.\n", + " hist, bins_x, bins_y = calgs.histogram2d(\n", + " mean_data.flatten().astype(np.float32),\n", + " pulseId.flatten().astype(np.float32),\n", + " bins=(100, int(pulse_range[1])),\n", + " range=[[clamp(vmin, -50, -0.2), min(vmax, 20000)], pulse_range]\n", + " )\n", + " do_2d_plot(hist, (bins_x, bins_y), \"Signal (ADU)\", \"Pulse id\", title=\"Signal-Pulse ID (Extended View)\")" ] }, { @@ -1143,7 +1150,7 @@ "fig = plt.figure(figsize=(10, 10))\n", "corrected_ave = np.nansum(corrected, axis=(2, 3))\n", "plt.scatter(corrected_ave.flatten()/10**6, blshift.flatten(), s=0.9)\n", - "plt.xlim(-1, 1000)\n", + "plt.xlim(np.nanpercentile(corrected_ave/10**6, [2, 98]))\n", "plt.grid()\n", "plt.xlabel('Illuminated corrected [MADU] ')\n", "_ = plt.ylabel('Estimated baseline shift [ADU]')" @@ -1176,8 +1183,9 @@ " fig = plt.figure(figsize=(20, 10))\n", " ax = fig.add_subplot(111)\n", " data = np.mean(raw[slice(*cell_sel.crange), 0, ...], axis=0)\n", - " vmin, vmax = get_range(data, 5)\n", - " ax = geom.plot_data_fast(data, ax=ax, cmap=\"jet\", vmin=vmin, vmax=vmax)\n", + " vmin, vmax = np.percentile(data, [5, 95])\n", + " ax = geom.plot_data_fast(data, ax=ax, vmin=vmin, vmax=vmax, cmap=cmap)\n", + " pass\n", "else:\n", " print(\"Skipping mean RAW preview for single memory cell, \"\n", " f\"see single shot image for selected cell ID {cell_id_preview}.\")" @@ -1192,8 +1200,10 @@ "display(Markdown(f'Single shot of the RAW data from cell {cell_id_preview} \\n'))\n", "fig = plt.figure(figsize=(20, 10))\n", "ax = fig.add_subplot(111)\n", - "vmin, vmax = get_range(raw[cell_idx_preview, 0, ...], 5)\n", - "ax = geom.plot_data_fast(raw[cell_idx_preview, 0, ...], ax=ax, cmap=\"jet\", vmin=vmin, vmax=vmax)" + "vmin, vmax = np.percentile(raw[cell_idx_preview, 0, ...], [5, 95])\n", + "ax = geom.plot_data_fast(\n", + " raw[cell_idx_preview, 0, ...], ax=ax, vmin=vmin, vmax=vmax, cmap=cmap)\n", + "pass" ] }, { @@ -1209,8 +1219,9 @@ " fig = plt.figure(figsize=(20, 10))\n", " ax = fig.add_subplot(111)\n", " data = np.mean(corrected, axis=0)\n", - " vmin, vmax = get_range(data, 7)\n", - " ax = geom.plot_data_fast(data, ax=ax, cmap=\"jet\", vmin=-50, vmax=vmax)\n", + " vmax = np.nanpercentile(data, 99.8)\n", + " ax = geom.plot_data_fast(data, ax=ax, vmin=max(-50, np.nanmin(data)), vmax=vmax, cmap=cmap)\n", + " pass\n", "else:\n", " print(\"Skipping mean CORRECTED preview for single memory cell, \"\n", " f\"see single shot image for selected cell ID {cell_id_preview}.\")" @@ -1225,9 +1236,15 @@ "display(Markdown(f'A single shot of the CORRECTED image from cell {cell_id_preview} \\n'))\n", "fig = plt.figure(figsize=(20, 10))\n", "ax = fig.add_subplot(111)\n", - "vmin, vmax = get_range(corrected[cell_idx_preview], 7, -50)\n", - "vmin = - 50\n", - "ax = geom.plot_data_fast(corrected[cell_idx_preview], ax=ax, cmap=\"jet\", vmin=vmin, vmax=vmax)" + "vmax = np.nanpercentile(corrected[cell_idx_preview], 99.8)\n", + "ax = geom.plot_data_fast(\n", + " corrected[cell_idx_preview],\n", + " ax=ax,\n", + " vmin=max(-50, np.nanmin(corrected[cell_idx_preview])),\n", + " vmax=vmax,\n", + " cmap=cmap,\n", + ")\n", + "pass" ] }, { @@ -1239,13 +1256,15 @@ "fig = plt.figure(figsize=(20, 10))\n", "ax = fig.add_subplot(111)\n", "vmin, vmax = get_range(corrected[cell_idx_preview], 5, -50)\n", - "nbins = np.int((vmax + 50) / 2)\n", + "nbins = int((vmax + 50) / 2)\n", "h = ax.hist(corrected[cell_idx_preview].flatten(),\n", " bins=nbins, range=(-50, vmax),\n", " histtype='stepfilled', log=True)\n", "plt.xlabel('[ADU]')\n", "plt.ylabel('Counts')\n", - "ax.grid()" + "ax.grid()\n", + "plt.title(f'Log-scaled histogram for corrected data for cell {cell_idx_preview}')\n", + "pass" ] }, { @@ -1259,8 +1278,8 @@ "vmin, vmax = get_range(corrected, 10, -100)\n", "vmax = np.nanmax(corrected)\n", "if vmax > 50000:\n", - " vmax=50000\n", - "nbins = np.int((vmax + 100) / 5)\n", + " vmax = 50000\n", + "nbins = int((vmax + 100) / 5)\n", "h = ax.hist(corrected.flatten(), bins=nbins,\n", " range=(-100, vmax), histtype='step', log=True, label = 'All')\n", "ax.hist(corrected[gains == 0].flatten(), bins=nbins, range=(-100, vmax),\n", @@ -1272,7 +1291,9 @@ "ax.legend()\n", "ax.grid()\n", "plt.xlabel('[ADU]')\n", - "plt.ylabel('Counts')" + "plt.ylabel('Counts')\n", + "plt.title(f'Overlaid Histograms for corrected data for multiple gains')\n", + "pass" ] }, { @@ -1293,8 +1314,10 @@ "source": [ "fig = plt.figure(figsize=(20, 10))\n", "ax = fig.add_subplot(111)\n", - "ax = geom.plot_data_fast(np.max(gains, axis=0), ax=ax,\n", - " cmap=\"jet\", vmin=-1, vmax=3)" + "ax = geom.plot_data_fast(\n", + " np.max(gains, axis=0), ax=ax,\n", + " cmap=cmap, vmin=-0.3, vmax=2.3) # Extend cmap for wrong gain values.\n", + "pass" ] }, { @@ -1336,7 +1359,9 @@ "source": [ "fig = plt.figure(figsize=(20, 10))\n", "ax = fig.add_subplot(111)\n", - "geom.plot_data_fast(np.log2(mask[cell_idx_preview]), ax=ax, vmin=0, vmax=32, cmap=\"jet\")" + "geom.plot_data_fast(\n", + " np.log2(mask[cell_idx_preview]), ax=ax, vmin=0, vmax=32, cmap=cmap)\n", + "pass" ] }, { @@ -1384,7 +1409,9 @@ "source": [ "fig = plt.figure(figsize=(20, 10))\n", "ax = fig.add_subplot(111)\n", - "geom.plot_data_fast(np.mean(mask>0, axis=0), vmin=0, ax=ax, vmax=1, cmap=\"jet\")" + "geom.plot_data_fast(\n", + " np.mean(mask>0, axis=0), vmin=0, ax=ax, vmax=1, cmap=cmap)\n", + "pass" ] }, { @@ -1404,8 +1431,9 @@ "ax = fig.add_subplot(111)\n", "cm = np.copy(mask)\n", "cm[cm > BadPixels.NO_DARK_DATA.value] = 0\n", - "ax = geom.plot_data_fast(np.mean(cm>0, axis=0),\n", - " vmin=0, ax=ax, vmax=1, cmap=\"jet\")" + "ax = geom.plot_data_fast(\n", + " np.mean(cm>0, axis=0), vmin=0, ax=ax, vmax=1, cmap=cmap)\n", + "pass" ] } ],