{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# LPD Radial X-ray Gain Evaluation #\n",
"\n",
"Author: S. Hauf, Version 0.5\n",
"\n",
"Taking proper flat field for LPD can be difficult, as air scattering will always be present. Additionally, the large detector mandates a large distance to the source, in order to reduce $1/r$ effects. \n",
"\n",
"Because of this a radial evaluation method is used, which assumes that pixels a the same radial distance $r$ should on average have the same signal $S(r)$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"in_folder = \"/gpfs/exfel/exp/FXE/201831/p900038/proc/\" # path to already corrected input data, required\n",
"out_folder = \"/gpfs/exfel/exp/FXE/201831/p900038/usr/calibration0818/FF/\" # path to output to, required\n",
"run = 125 # runs to process, required\n",
"sequences = [0] # which sequence files to use\n",
"capacitor_setting = 5 # capacitor_setting for which data was taken\n",
"\n",
"mem_cells = 512 # number of memory cells used\n",
"local_output = True # output constants locally\n",
"db_output = False # output constants to database\n",
"bias_voltage = 250 # detector bias voltage\n",
"cal_db_interface = \"tcp://max-exfl015:5005\" # the database interface to use\n",
"\n",
"use_dir_creation_date = True # use the creation date of the directory for database time derivation\n",
"instrument = \"FXE\"\n",
"limit_trains = 10 # limit the number of train for the evaluation\n",
"geometry_file = \"/gpfs/exfel/d/cal/exchange/lpdMF_00.h5\" # the geometry file to use, MAR 2018\n",
"beam_center_offset = [1.5, 1] # offset from the beam center, MAR 2018\n",
"allowed_gain_thresholds = [0.5, 1.5] # relative gain values within these bounds are valid\n",
"badpix_entire_asic_threshold = 0.25 # if more than this fraction of pixels on an ASIC are \"bad\" the entire ASIC is flagged\n",
"photon_energy = 9.2 # photon enery in keV\n",
"max_images = 1024 # maximum number of images to use in evaluation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import warnings\n",
"\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"# make sure a cluster is running with ipcluster start --n=32, give it a while to start\n",
"\n",
"import os\n",
"\n",
"import h5py\n",
"import matplotlib\n",
"import numpy as np\n",
"\n",
"matplotlib.use(\"agg\")\n",
"import matplotlib.pyplot as plt\n",
"\n",
"%matplotlib inline\n",
"\n",
"\n",
"from collections import OrderedDict\n",
"from datetime import datetime\n",
"\n",
"import cal_tools.metrology as metro\n",
"from cal_tools.enums import BadPixels\n",
"from cal_tools.plotting import create_constant_overview, plot_badpix_3d, show_overview\n",
"from cal_tools.tools import (\n",
" gain_map_files,\n",
" get_constant_from_db,\n",
" get_dir_creation_date,\n",
" get_notebook_name,\n",
" parse_runs,\n",
" run_prop_seq_from_path,\n",
")\n",
"from iCalibrationDB import Conditions, ConstantMetaData, Constants, Detectors, Versions\n",
"\n",
"creation_time = None\n",
"if use_dir_creation_date:\n",
" creation_time = get_dir_creation_date(in_folder, run)\n",
" print(\"Using {} as creation time\".format(creation_time))\n",
"\n",
"run = \"r{:04d}\".format(run)\n",
"in_folder = \"{}/{}\".format(in_folder, run)\n",
"\n",
"gains = np.arange(3)\n",
"max_cells = mem_cells\n",
"cells = np.arange(max_cells)\n",
"\n",
"\n",
"QUADRANTS = 4\n",
"MODULES_PER_QUAD = 4\n",
"DET_FILE_INSET = \"LPD\"\n",
"\n",
"cap = \"{}pf\".format(capacitor_setting)\n",
"\n",
"\n",
"if not os.path.exists(out_folder):\n",
" os.makedirs(out_folder)\n",
"\n",
"run, proposal, seq = run_prop_seq_from_path(in_folder)\n",
"\n",
"# start lower right and then counter-clock-wise\n",
"dc = beam_center_offset\n",
"\n",
"#d_quads = [(-15,-300),(10,-10),(-255,15),(-280,-275)]\n",
"#d_quads = [(10,-275),(-15,15),(-280,-10),(-255,-300)] # DEC 2017\n",
"d_quads = [(-14+dc[0],-300+dc[1]),(10+dc[0],-9+dc[1]),(-256+dc[0],15+dc[1]),(-280+dc[0],-276+dc[1])] # MAR 2018\n",
"#d_quads = [(-16+dc[0],-300+dc[1]),(10+dc[0],-9+dc[1]),(-256+dc[0],15+dc[1]),(-280+dc[0],-276+dc[1])] # MAY 2018"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# set everything up filewise\n",
"from queue import Queue\n",
"\n",
"if not os.path.exists(out_folder):\n",
" os.makedirs(out_folder)\n",
" \n",
"def map_modules_from_files(filelist):\n",
" module_files = {}\n",
" mod_ids = {}\n",
" for quadrant in range(0, QUADRANTS):\n",
" for module in range(0, MODULES_PER_QUAD):\n",
" name = \"Q{}M{}\".format(quadrant + 1, module + 1)\n",
" module_files[name] = Queue()\n",
" num = quadrant * 4 + module\n",
" mod_ids[name] = num\n",
" file_infix = \"{}{:02d}\".format(DET_FILE_INSET, num)\n",
" for file in filelist:\n",
" if file_infix in file:\n",
" module_files[name].put(file)\n",
" \n",
" return module_files, mod_ids\n",
"\n",
"dirlist = os.listdir(in_folder)\n",
"file_list = []\n",
"for entry in dirlist:\n",
" #only h5 file\n",
" abs_entry = \"{}/{}\".format(in_folder, entry)\n",
" if os.path.isfile(abs_entry) and os.path.splitext(abs_entry)[1] == \".h5\":\n",
" \n",
" if sequences is None:\n",
" file_list.append(abs_entry)\n",
" else:\n",
" for seq in sequences:\n",
" if \"{:05d}.h5\".format(seq) in abs_entry:\n",
" file_list.append(os.path.abspath(abs_entry))\n",
" \n",
"mapped_files, mod_ids = map_modules_from_files(file_list)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import copy\n",
"\n",
"import tabulate\n",
"from IPython.display import HTML, Latex, Markdown, display\n",
"\n",
"table = []\n",
"mfc = copy.copy(mapped_files)\n",
"ti = 0\n",
"for k, files in mfc.items():\n",
" i = 0\n",
" while not files.empty():\n",
" f = files.get()\n",
" if i == 0:\n",
" table.append((ti, k, i, f))\n",
" else:\n",
" table.append((ti, \"\", i, f))\n",
" i += 1\n",
" ti += 1\n",
"md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=[\"#\", \"module\", \"# module\", \"file\"]))) \n",
"mapped_files, mod_ids = map_modules_from_files(file_list)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"first_files = []\n",
"for i in range(16):\n",
" qm = \"Q{}M{}\".format(i//4 +1, i % 4 + 1)\n",
" if qm in mapped_files and not mapped_files[qm].empty():\n",
" fname_in = str(mapped_files[qm].get())\n",
" else:\n",
" print(\"Skipping {}\".format(qm))\n",
" first_files.append((None, None))\n",
" continue\n",
" fout = os.path.abspath(\"{}/{}\".format(in_folder, (os.path.split(fname_in)[-1]).replace(\"RAW\", \"CORR\")))\n",
" \n",
" first_files.append((fname_in, fout))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the flat fields all correction up to relative gain are applied. If for a given module no data exists, zero-length data is used instead."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"corrected = []\n",
"raw = []\n",
"gains = []\n",
"\n",
"for i, ff in enumerate(first_files):\n",
" try:\n",
" rf, cf = ff\n",
" \n",
" if rf is not None:\n",
" \n",
" infile = h5py.File(rf, \"r\")\n",
" raw.append(np.array(infile[\"/INSTRUMENT/FXE_DET_LPD1M-1/DET/{}CH0:xtdf/image/data\".format(i)][:max_images,0,...]))\n",
" infile.close()\n",
"\n",
" infile = h5py.File(cf, \"r\")\n",
" corrected.append(np.array(infile[\"/INSTRUMENT/FXE_DET_LPD1M-1/DET/{}CH0:xtdf/image/data\".format(i)][:max_images,...]))\n",
" gains.append(np.array(infile[\"/INSTRUMENT/FXE_DET_LPD1M-1/DET/{}CH0:xtdf/image/gain\".format(i)][:max_images,...]))\n",
" infile.close()\n",
" else:\n",
" raise\n",
" \n",
" except Exception as e:\n",
" corrected.append(np.zeros((max_images, 256, 256)))\n",
" raw.append(np.zeros((max_images, 256, 256)))\n",
" gains.append(np.zeros((max_images, 256, 256)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data Inspection ##\n",
"\n",
"It is important that the intensity in the image is radially symetric, and shows now strong intensity gradiants, e.g. through pronounced diffraction rings."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import cal_tools.metrology as metro\n",
"\n",
"in_files = \"{}/CORR*S{:05d}*.h5\".format(in_folder, sequences[0] if sequences else 0)\n",
"datapath = \"INSTRUMENT/FXE_DET_LPD1M-1/DET/{}CH0:xtdf/image/data\"\n",
"posarr = metro.positionFileList(in_files, datapath, geometry_file, d_quads, nImages = 100)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(15,15))\n",
"ax = fig.add_subplot(111)\n",
"posarr[posarr > 1e6] = np.nan\n",
"posarr[posarr < -1e4] = np.nan\n",
"parr = np.nanmean(posarr, axis=0)\n",
"im=ax.imshow((parr), vmin=0, vmax=max(3*np.median(parr[parr > 0]), 100))\n",
"cb = fig.colorbar(im)\n",
"cb.set_label(\"Intensity (ADU)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Additionally, data should ideally only be present in the highest gain stage. This can be verfied in the following plot, showing the maximum pixel gain value for the first 100 images."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"datapath = \"INSTRUMENT/FXE_DET_LPD1M-1/DET/{}CH0:xtdf/image/gain\"\n",
"posarr = metro.positionFileList(in_files, datapath, geometry_file, d_quads, nImages = 100)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(15,15))\n",
"ax = fig.add_subplot(111)\n",
"parr = np.nanmax(posarr, axis=0)\n",
"im=ax.imshow((parr), vmin=0, vmax=3)\n",
"cb = fig.colorbar(im)\n",
"cb.set_label(\"Gain Bit\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# create median value maps\n",
"means = []\n",
"max_gains = []\n",
"for i, c in enumerate(corrected):\n",
" #c[c>2**12] = np.nan\n",
" means.append(np.nanmean(c, axis=0))\n",
" max_gains.append(np.nanmax(gains[i], axis=0))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def splitChannelDataIntoTiles(channelData, clockwiseOrder=False):\n",
" \"\"\"Splits the raw channel data into indiviual tiles\n",
" \n",
" Args\n",
" ----\n",
" \n",
" channelData : ndarray\n",
" Raw channel data. Must have shape (256, 256)\n",
" \n",
" clockwiseOrder : bool, optional\n",
" If set to True, the sequence of tiles is given\n",
" in the clockwise order starting with the top\n",
" right tile (LPD standard). If set to false, tile\n",
" data is returned in reading order\n",
" \n",
" Returns\n",
" -------\n",
" \n",
" ndarray\n",
" Same data, but reshaped into (12, 32, 128)\n",
" \"\"\"\n",
" tiles = np.asarray(np.split(channelData, 8, axis=1))\n",
" tiles = np.asarray(np.split(tiles, 2, axis=1))\n",
" orderedTiles = np.moveaxis(tiles.reshape(16, 128, 32), 2, 1)\n",
" if clockwiseOrder:\n",
" # Naturally, the tile data after splitting is in reading\n",
" # order (i.e. top left tile is first, top right tile is second, \n",
" # etc.). The official LPD tile order however is clockwise, \n",
" # starting with the top right tile. The following array\n",
" # contains indices of tiles in reading order as they would\n",
" # be iterated in clockwise order (starting from the top right)\n",
" readingOrderToClockwise = list(range(8,16))+list(range(7,-1,-1))\n",
" # Return tiles in reading order\n",
" orderedTiles = orderedTiles[readingOrderToClockwise]\n",
" return orderedTiles"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"positions = []\n",
"mn_tiles = []\n",
"gn_tiles = []\n",
"\n",
"tile_order = [1,2,3,4]\n",
"for sm, mn in enumerate(means):\n",
" position = np.asarray([metro.getModulePosition(geometry_file,\n",
" 'Q4/M{:d}/T{:02d}'.format(tile_order[sm%4], idx+1))\n",
" for idx in range(16)])\n",
" positions.append(position)\n",
" mn_tile = metro.splitChannelDataIntoTiles(mn[::-1,::-1], clockwiseOrder=True)\n",
" mn_tiles.append(mn_tile)\n",
" \n",
" gn_tile = metro.splitChannelDataIntoTiles(max_gains[sm][::-1,::-1], clockwiseOrder=True)\n",
" gn_tiles.append(gn_tile)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from matplotlib.colors import LogNorm, PowerNorm\n",
"\n",
"\n",
"def translateToModuleBL(tilePositions):\n",
" tileHeight = 17.7 # mm\n",
" # The module origin is the top left corner of the\n",
" # top left tile. So in the clockwise order of LPD\n",
" # tiles, it is the last tile in the list\n",
" moduleOrigin = tilePositions[-1]\n",
" # np.asarray([1., -1.]): inverts the y-axis, since\n",
" # matplotlib.figure coordinate system is on the\n",
" # bottom left pointing up and to the right\n",
" moduleCoords = np.asarray([1., -1.])*(\n",
" # np.asarray([0., tileHeight]): Tile positions\n",
" # are translated from the top left corner to the\n",
" # bottom left corner, i.e. 0 in x, tileHeight in y\n",
" tilePositions + np.asarray([0., tileHeight]) - moduleOrigin\n",
" )\n",
" # In the clockwise order of LPD tiles, the 8th\n",
" # tile in the list is the bottom left tile\n",
" bottomLeft8th = np.asarray([0., moduleCoords[8][1]])\n",
" # Translate coordinates to the bottom left corner \n",
" # of the bottom left tile\n",
" bottomLeft = moduleCoords - bottomLeft8th\n",
" return bottomLeft\n",
"\n",
"\n",
"def plotSupermoduleData(smData, smPositions, dquads, zoom=1., vmin=1., vmax=2**16*0.5, norm=\"power\", gamma=0.3):\n",
" \n",
" # Conversion coefficient, required since\n",
" # matplotlib does its business in inches\n",
" mmToInch = 1./25.4 # inch/mm\n",
" \n",
" # Some constants\n",
" \n",
" numberOfRows = 16\n",
" numberOfCols = 4\n",
" tileWidth = 65.7 # in mm\n",
" tileHeight = 17.7 # in mm\n",
" \n",
" # Base width and height are given by spatial\n",
" # extend of the modules. The constants 3.4 and 1\n",
" # are estimated as a best guess for gaps between\n",
" # modules.\n",
" figureWidth = zoom * numberOfCols*(tileWidth + 3.4)*mmToInch\n",
" figureHeight = zoom * numberOfRows*(tileHeight + 1.)*mmToInch\n",
" fig = plt.figure(figsize=(figureWidth, figureHeight))\n",
" \n",
" quads = []\n",
" quad_pos = []\n",
" for q in range(len(smData)//4):\n",
" dq = min(4, len(smData)-4*q)\n",
" qdata = np.concatenate(smData[4*q:4*q+dq], axis=0)\n",
" qpos = np.concatenate(smPositions[4*q:4*q+dq], axis=0)\n",
" quads.append(qdata)\n",
" quad_pos.append(qpos)\n",
" \n",
" \n",
" for q, tileData in enumerate(quads):\n",
" \n",
" metrologyPositions = quad_pos[q]\n",
" numberOfTiles = tileData.shape[0]\n",
" \n",
" # The metrology file references module positions \n",
" bottomRightCornerCoordinates = metrologyPositions #translateToModuleBL(metrologyPositions)\n",
" \n",
" # The offset here accounts for the fact that there\n",
" # might be negative x,y values\n",
" offset = np.asarray(\n",
" [min(bottomRightCornerCoordinates[:, 0]), \n",
" min(bottomRightCornerCoordinates[:, 1])]\n",
" )\n",
"\n",
" # Account for blank borders in the plot\n",
" borderLeft = 0.5 * mmToInch\n",
" borderBottom = 0.5 * mmToInch\n",
"\n",
" # The height and width of the plot remain\n",
" # constant for a given supermodule\n",
" width = zoom * 65.7 * mmToInch / (figureWidth - 2.*borderLeft)\n",
" height = zoom * 17.7 * mmToInch / (figureHeight - 2.*borderBottom)\n",
"\n",
" for i in range(numberOfTiles):\n",
" # This is the top left corner of the tile with\n",
" # respect to the top left corner of the supermodule\n",
" x0, y0 = bottomRightCornerCoordinates[i] + dquads[q]\n",
" # Transform to figure coordinates\n",
" ax0 = borderLeft + zoom * x0 * mmToInch / (figureWidth - 2.*borderLeft)\n",
" ay0 = borderBottom + zoom * y0 * mmToInch / (figureHeight - 2.*borderBottom)\n",
" # Place the plot in the figure\n",
" ax = fig.add_axes((ax0, ay0, width, height), frameon=False)\n",
" # Do not display axes, tick markers or labels\n",
" ax.tick_params(\n",
" axis='both', left='off', top='off', right='off', bottom='off', \n",
" labelleft='off', labeltop='off', labelright='off', labelbottom='off'\n",
" )\n",
" # Plot the image. If one wanted to have a colorbar\n",
" # the img object would be needed to produce one\n",
" if norm ==\"power\":\n",
" img = ax.imshow(\n",
" tileData[i][::-1,:], \n",
" interpolation='nearest', \n",
" norm=PowerNorm(vmin=vmin, vmax=vmax, gamma=gamma)\n",
" )\n",
" else:\n",
" img = ax.imshow(\n",
" tileData[i][::-1,:], \n",
" interpolation='nearest', \n",
" vmin=vmin, vmax=vmax\n",
" )\n",
" #fig.subplots_adjust(right=0.8)\n",
" cbar_ax = fig.add_axes([1.05, 0.15, 0.05, 0.7])\n",
" cb = fig.colorbar(img, cax=cbar_ax)\n",
" cb.ax.tick_params(labelsize=40) \n",
" \n",
" \n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# convert the Carthesian coordinates of the detector to polar coordinates\n",
"def mod_cart_to_pol(d, dx, dy, filter_by_val=True):\n",
" \"\"\" Convert carthesian coords to polar coords\n",
" \"\"\"\n",
" cx, cy = d.shape\n",
" x = np.arange(cx)+dx\n",
" y = np.arange(cy)+dy\n",
" x = np.repeat(x[:,None], cy, axis=1)\n",
" y = np.repeat(y[None,:], cx, axis=0)\n",
" \n",
" rho = np.sqrt(x**2 + y**2).flatten()\n",
" phi = np.arctan2(y, x).flatten()\n",
" flat = d.flatten()\n",
" \n",
" # we also perform a bit of filtering here\n",
" \n",
" if filter_by_val:\n",
" good = np.isfinite(flat) & (flat > 0) & (flat < 1e5)\n",
" \n",
" return rho[good], phi[good], flat[good], good\n",
" \n",
" return rho, phi, flat, None"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Here we create histograms of the data in a polar coordinate system.\n",
"# We use numpys hist2d function, giving it the polar coordinates of \n",
"# each pixels, and weighing that coordinate with the pixel's value.\n",
"# We obtain a histogram for each module, according to its position defined\n",
"# in the coord_list.\n",
"from scipy.stats import binned_statistic_2d\n",
"\n",
"hs = []\n",
"bins_nums = []\n",
"edges = []\n",
"#d_quads = [(-15,-300),(10,-10),(-255,15),(-280,-275)]\n",
"#d_quads = [(15,-275),(-15,15),(-280,-10),(-250,-300)]\n",
"goods = []\n",
"bins = 5000\n",
"for q in range(4):\n",
" quads = []\n",
" quad_pos = []\n",
"\n",
" dq = 4\n",
" qdata = np.concatenate(mn_tiles[4*q:4*q+dq], axis=0)\n",
" qpos = np.concatenate(positions[4*q:4*q+dq], axis=0)\n",
" \n",
" for s in range(qdata.shape[0]):\n",
" \n",
" \n",
" dx, dy = qpos[s,...].tolist()\n",
" dx += d_quads[q][0]\n",
" dy += d_quads[q][1]\n",
" dy /= 0.5\n",
" dx /= 0.5\n",
"\n",
" td = qdata[s,...]\n",
" rho, phi, weights, good = mod_cart_to_pol(td, dy, dx)\n",
" #h, phi_edges, rho_edges = np.histogram2d(phi, rho, bins=(1000,1000),\n",
" # range=((-np.pi, np.pi), (0, 1000)),\n",
" # weights=weights)\n",
" h, phi_edges, rho_edges, bns = binned_statistic_2d(phi, rho, weights, bins=(bins,bins),\n",
" range=((-np.pi, np.pi), (0, 1000)),\n",
" statistic = \"sum\")\n",
" bins_nums.append(bns)\n",
" hs.append(h)\n",
" edges.append((phi_edges, rho_edges))\n",
" goods.append(good)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Polar Plot ##\n",
"\n",
"Plotting the data in polar coordinates should show an intensity gradiant from left to right with regions of similar intensity extending vertically."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Plot data in polar coordinates. Diffraction rings should now show as vertical\n",
"# lines, connected over all modules - if the positions in the coordinate list\n",
"# are indeed correct that is.\n",
"fig = plt.figure(figsize=(20,20))\n",
"ax = fig.add_subplot(111)\n",
"x = np.arange(bins)/bins*1000*500e-6\n",
"y = np.arange(bins)/bins*2.\n",
"ds = np.array(hs).sum(axis=0)\n",
"im = plt.pcolormesh(x, y, ds, vmin=0, vmax=5e2)\n",
"cb = fig.colorbar(im)\n",
"ax.set_xlim(0, np.max(x))\n",
"ax.set_ylim(np.min(y), np.max(y))\n",
"ax.set_xlabel(\"r (m)\")\n",
"ax.set_ylabel(\"$\\phi$ ($\\pi$)\")\n",
"cb.set_label(\"Intensity (arb. units)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Accordingly, an azimutally integrated profile should should be smooth, without any sharp peaks or artifacts. Additionally, the mean intensity should be at least 150 ADU in all reagions."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [],
"source": [
"# With appropriate coordinates given, plotting a profile along the\n",
"# vertical axis should give us the positions of the diffraction peaks,\n",
"# Here still as distances on the detector plane. With knowledge of the\n",
"# detector to sample distance, these could then be converted in \n",
"# reciprocal coordinates.\n",
"ds[ds == 0] = np.nan\n",
"profile = np.nanmedian(ds, axis=0)\n",
"fig = plt.figure(figsize=(15,5))\n",
"ax = fig.add_subplot(111)\n",
"ax.plot(x, profile)\n",
"ax.set_xlabel(\"r (m)\")\n",
"ax.set_ylabel(\"Median intensity (arb. units)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Relative Gain Calculation ##\n",
"\n",
"The relative gain is then given by each (polar) pixels intensity w.r.t to the mean instensity of all pixels at the same radial distance."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"relgain = ds/profile[None,:]\n",
"\n",
"fig = plt.figure(figsize=(20,20))\n",
"ax = fig.add_subplot(111)\n",
"x = np.arange(bins)/bins*1000*500e-6\n",
"y = np.arange(bins)/bins*2.\n",
"ds = np.array(hs).sum(axis=0)\n",
"im = plt.pcolormesh(x, y, relgain, vmin=0.8, vmax=1.2)\n",
"cb = fig.colorbar(im)\n",
"ax.set_xlim(0, np.max(x))\n",
"ax.set_ylim(np.min(y), np.max(y))\n",
"ax.set_xlabel(\"r (m)\")\n",
"ax.set_ylabel(\"$\\phi$ ($\\pi$)\")\n",
"cb.set_label(\"Intensity (arb. units)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Correcting the input data witht this gain map should yield a homogenoeous image in polar coordinates"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20,20))\n",
"ax = fig.add_subplot(111)\n",
"x = np.arange(bins)/bins*1000*500e-6\n",
"y = np.arange(bins)/bins*2.\n",
"ds = np.array(hs).sum(axis=0)\n",
"im = plt.pcolormesh(x, y, ds/relgain, vmin=0, vmax=1e3)\n",
"cb = fig.colorbar(im)\n",
"ax.set_xlim(0, np.max(x))\n",
"ax.set_ylim(np.min(y), np.max(y))\n",
"ax.set_xlabel(\"r (m)\")\n",
"ax.set_ylabel(\"$\\phi$ ($\\pi$)\")\n",
"cb.set_label(\"Intensity (arb. units)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, the gain map can be tranferred back into carthesian coordinates, and re-attributed to the individual tiles and supermodules"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"i = 0\n",
"rgain_tiles = []\n",
"for q in range(4):\n",
" quads = []\n",
" quad_pos = []\n",
"\n",
" dq = 4\n",
" qdata = np.zeros_like(np.concatenate(mn_tiles[4*q:4*q+dq], axis=0)) \n",
" qpos = np.concatenate(positions[4*q:4*q+dq], axis=0)\n",
" \n",
" for s in range(qdata.shape[0]):\n",
" \n",
" dx, dy = qpos[s,...].tolist()\n",
" dx += d_quads[q][0]\n",
" dy += d_quads[q][1]\n",
" dy /= 0.5\n",
" dx /= 0.5\n",
" \n",
" rgain = hs[i]/profile[None,:]\n",
" good = goods[i]\n",
" this_tile = np.zeros_like(qdata[s,...])\n",
" \n",
" rho, phi, _, _ = mod_cart_to_pol(this_tile, dy, dx, filter_by_val=False)\n",
" # bin edges will contain lower and upper limits, which go past the range values\n",
" x_inds = np.searchsorted(phi_edges[1:-1], phi[good])\n",
" y_inds = np.searchsorted(rho_edges[1:-1] , rho[good])\n",
" #x_inds[x_inds >= rgain.shape[0]] = rgain.shape[0]-1\n",
" #y_inds[y_inds >= rgain.shape[1]] = rgain.shape[1]-1\n",
" \n",
" this_tile.flat[good] = rgain[x_inds, y_inds].flat\n",
" \n",
" qdata[s,...] = this_tile\n",
" i += 1\n",
" \n",
" rgain_tiles += np.split(qdata, 4, axis=0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plotSupermoduleData(rgain_tiles, positions, d_quads, vmin=0.5, vmax=1.5, norm=\"lin\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"corr_tiles = []\n",
"photon_gain = 22.4 # ADU per photon\n",
"for i, tile in enumerate(mn_tiles):\n",
" ctile = tile/rgain_tiles[i]/photon_gain\n",
" corr_tiles.append(ctile)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"max_corr_tiles = []\n",
"mean_corr_tiles = []\n",
"for m, sm in enumerate(corrected):\n",
" ind_cor = []\n",
" for i in range(sm.shape[0]):\n",
" tiled = metro.splitChannelDataIntoTiles(sm[i, ::-1,::-1], clockwiseOrder=True)\n",
" tiled /= rgain_tiles[m]\n",
" ind_cor.append(tiled/photon_gain)\n",
" max_corr_tiles.append(np.nanmax(np.array(ind_cor), axis=0))\n",
" mean_corr_tiles.append(np.nanmean(np.array(ind_cor), axis=0))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Verification ##\n",
"\n",
"Plotting the maximum pixel value of the input data, corrected with the gain map derived from it will enhance any artifacts."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plotSupermoduleData(max_corr_tiles, positions, d_quads, vmin=0, vmax=5*np.nanmedian(max_corr_tiles), gamma=0.7)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mean pixel value ###"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plotSupermoduleData(mean_corr_tiles, positions, d_quads, vmin=0, vmax=1.5*np.nanmedian(max_corr_tiles), gamma=0.7)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# now arrange back to supermodules\n",
"sm_gains = []\n",
"for i, rgain in enumerate(rgain_tiles):\n",
" sm_gain = np.zeros((256,256))\n",
" for j in range(8):\n",
" sm_gain[j*32:(j+1)*32,128:] = rgain[j,...]\n",
" sm_gain[256-(j+1)*32:256-j*32,:128] = rgain[j+8,...]\n",
" sm_gains.append(sm_gain)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bad Pixels ##\n",
"\n",
"Bad pixels are defined by thresholds and non-converging evaluation. Additionally, full ASICs are masked if a significant fraction of pixels on a given ASIC is considered \"bad\"."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"bad_pixels = []\n",
"num_asics_discarded = 0\n",
"full_modules_missing = 0\n",
"bpix_total = 0\n",
"for i, gain in enumerate(sm_gains):\n",
" bpix = np.zeros(gain.shape, np.uint32)\n",
" bpix[~np.isfinite(gain)] |= BadPixels.FF_GAIN_EVAL_ERROR.value\n",
" bpix[gain == 0] |= BadPixels.FF_NO_ENTRIES.value\n",
" bpix[(gain < allowed_gain_thresholds[0]) | (gain > allowed_gain_thresholds[1])] |= BadPixels.FF_GAIN_DEVIATION.value\n",
" for j in range(8):\n",
" for k in range(16):\n",
" asic_data = bpix[j*32:(j+1)*32, k*16:(k+1)*16]\n",
" bp_frac = float(np.count_nonzero(asic_data))/asic_data.size\n",
" \n",
" if bp_frac > badpix_entire_asic_threshold:\n",
" num_asics_discarded += 1\n",
" bpix[j*32:(j+1)*32, k*16:(k+1)*16] |= BadPixels.FF_GAIN_DEVIATION.value\n",
" if np.all(bpix != 0):\n",
" full_modules_missing += 1\n",
" bpix_total += np.count_nonzero(bpix) \n",
" bad_pixels.append(bpix)\n",
"print(\"ASICs with bad pixel fraction over {}: {}\".format(badpix_entire_asic_threshold,\n",
" num_asics_discarded-full_modules_missing*16*8))\n",
"print(\"Full modules considered bad: {}\".format(full_modules_missing))\n",
"print(\"Bad pixel fraction (excluding missing modules): {:0.2f}\".format((bpix_total-full_modules_missing*256*256)/\n",
" (1024*1024-full_modules_missing*256*256)))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"bad_pixels_tiles = []\n",
"for m, sm in enumerate(bad_pixels):\n",
" bp = []\n",
" \n",
" tiled = metro.splitChannelDataIntoTiles(sm[::-1,::-1], clockwiseOrder=True)\n",
" bp.append(np.log2(tiled))\n",
" bad_pixels_tiles.append(np.nanmax(np.array(bp), axis=0))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plotSupermoduleData(bad_pixels_tiles, positions, d_quads, vmin=0, vmax=32)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ofile = \"{}/lpd_flatfield_store_{}.h5\".format(out_folder, run)\n",
"store_file = h5py.File(ofile, \"w\")\n",
"for i, sgain in enumerate(sm_gains):\n",
" qm = \"Q{}M{}\".format(i // 4 + 1, i % 4 +1)\n",
" store_file[\"{}/FlatField/0/data\".format(qm)] = sgain \n",
" store_file[\"{}/BadPixelsFF/0/data\".format(qm)] = bad_pixels[i] \n",
"store_file.close()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"if db_output:\n",
" for i, sgain in enumerate(sm_gains):\n",
" qm = \"Q{}M{}\".format(i // 4 + 1, i % 4 +1)\n",
" \n",
" metadata = ConstantMetaData() \n",
" slopes = Constants.LPD.SlopesFF()\n",
" slopes.data = sgain\n",
" metadata.calibration_constant = slopes\n",
"\n",
" # set the operating condition\n",
" condition = Conditions.Illuminated.LPD(1., bias_voltage, photon_energy,\n",
" capacitor=capacitor_setting, beam_energy=None)\n",
" device = getattr(Detectors.LPD1M1, qm)\n",
"\n",
" if device:\n",
"\n",
" metadata.detector_condition = condition\n",
"\n",
" # specify the a version for this constant\n",
" if creation_time is None:\n",
" metadata.calibration_constant_version = Versions.Now(device=device)\n",
" else:\n",
" metadata.calibration_constant_version = Versions.Timespan(device=device,\n",
" start=creation_time)\n",
" metadata.send(cal_db_interface)\n",
" \n",
" metadata = ConstantMetaData() \n",
" badpix = Constants.LPD.BadPixelsFF()\n",
" badpix.data = bad_pixels[i] \n",
" metadata.calibration_constant = badpix\n",
"\n",
" # set the operating condition\n",
" condition = Conditions.Illuminated.LPD(1., bias_voltage, photon_energy,\n",
" capacitor=capacitor_setting, beam_energy=None)\n",
" device = getattr(Detectors.LPD1M1, qm)\n",
"\n",
" if device:\n",
" metadata.detector_condition = condition\n",
"\n",
" # specify the a version for this constant\n",
" if creation_time is None:\n",
" metadata.calibration_constant_version = Versions.Now(device=device)\n",
" else:\n",
" metadata.calibration_constant_version = Versions.Timespan(device=device,\n",
" start=creation_time)\n",
" metadata.send(cal_db_interface)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"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.6"
}
},
"nbformat": 4,
"nbformat_minor": 1
}