{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# DSSC Offline Correction #\n",
"\n",
"Author: European XFEL Detector Group, Version: 1.0\n",
"\n",
"Offline Calibration for the DSSC Detector"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-21T11:30:06.730220Z",
"start_time": "2019-02-21T11:30:06.658286Z"
},
"collapsed": true
},
"outputs": [],
"source": [
"in_folder = \"/gpfs/exfel/exp/SCS/201802/p002252/raw/\" # the folder to read data from, required\n",
"run = 20 # runs to process, required\n",
"out_folder = \"/gpfs/exfel/data/scratch/xcal/test/\" # the folder to output to, required\n",
"sequences = [-1] # sequences to correct, set to -1 for all, range allowed\n",
"mem_cells = 0 # number of memory cells used, set to 0 to automatically infer\n",
"overwrite = True # set to True if existing data should be overwritten\n",
"cluster_profile = \"noDB\" # cluster profile to use\n",
"max_pulses = 500 # maximum number of pulses per train\n",
"bias_voltage = 300 # detector bias voltage\n",
"cal_db_interface = \"tcp://max-exfl016:8020#8025\" # the database interface to use\n",
"use_dir_creation_date = True # use the creation data of the input dir for database queries\n",
"sequences_per_node = 2 # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel\n",
"cal_db_timeout = 30000 # in milli seconds\n",
"chunk_size_idim = 1 # chunking size of imaging dimension, adjust if user software is sensitive to this.\n",
"instrument = \"SCS\" # the instrument the detector is installed at, required\n",
"mask_noisy_asic = 0.25 # set to a value other than 0 and below 1 to mask entire ADC if fraction of noisy pixels is above\n",
"offset_image = \"-1\" # last one\n",
"mask_cold_asic = 0.25 # mask cold ASICS if number of pixels with negligable standard deviation is larger than this fraction\n",
"noisy_pix_threshold = 1. # threshold above which ap pixel is considered noisy.\n",
"geo_file = \"/gpfs/exfel/data/scratch/xcal/dssc_geo_june19.h5\" # detector geometry file\n",
"\n",
"\n",
"def balance_sequences(in_folder, run, sequences, sequences_per_node):\n",
" import glob\n",
" import re\n",
" import numpy as np\n",
" if sequences[0] == -1:\n",
" sequence_files = glob.glob(\"{}/r{:04d}/*-S*.h5\".format(in_folder, run))\n",
" seq_nums = set()\n",
" for sf in sequence_files:\n",
" seqnum = re.findall(r\".*-S([0-9]*).h5\", sf)[0]\n",
" seq_nums.add(int(seqnum))\n",
" seq_nums -= set(sequences)\n",
" else:\n",
" seq_nums = set(sequences)\n",
" nsplits = len(seq_nums)//sequences_per_node+1\n",
" while nsplits > 32:\n",
" sequences_per_node += 1\n",
" nsplits = len(seq_nums)//sequences_per_node+1\n",
" print(\"Changed to {} sequences per node to have a maximum of 8 concurrent jobs\".format(sequences_per_node))\n",
" return [l.tolist() for l in np.array_split(list(seq_nums), nsplits) if l.size > 0]\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-21T11:30:07.086286Z",
"start_time": "2019-02-21T11:30:06.929722Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"import sys\n",
"from collections import OrderedDict\n",
"\n",
"# make sure a cluster is running with ipcluster start --n=32, give it a while to start\n",
"import os\n",
"import h5py\n",
"import numpy as np\n",
"import matplotlib\n",
"matplotlib.use(\"agg\")\n",
"import matplotlib.pyplot as plt\n",
"from ipyparallel import Client\n",
"print(\"Connecting to profile {}\".format(cluster_profile))\n",
"view = Client(profile=cluster_profile)[:]\n",
"view.use_dill()\n",
"\n",
"from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions\n",
"from cal_tools.tools import (gain_map_files, parse_runs, run_prop_seq_from_path, get_notebook_name,\n",
" get_dir_creation_date, get_constant_from_db)\n",
"\n",
"from dateutil import parser\n",
"from datetime import timedelta\n",
"\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",
"in_folder = \"{}/r{:04d}\".format(in_folder, run)\n",
"\n",
"\n",
"if sequences[0] == -1:\n",
" sequences = None\n",
" \n",
"\n",
"QUADRANTS = 4\n",
"MODULES_PER_QUAD = 4\n",
"DET_FILE_INSET = \"DSSC\"\n",
"CHUNK_SIZE = 512\n",
"MAX_PAR = 32\n",
"\n",
"if in_folder[-1] == \"/\":\n",
" in_folder = in_folder[:-1]\n",
"out_folder = \"{}/{}\".format(out_folder, os.path.split(in_folder)[-1])\n",
"print(\"Outputting to {}\".format(out_folder))\n",
"\n",
"if not os.path.exists(out_folder):\n",
" os.makedirs(out_folder)\n",
"elif not overwrite:\n",
" raise AttributeError(\"Output path exists! Exiting\")\n",
"\n",
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"\n",
"loc = None\n",
"if instrument == \"SCS\":\n",
" loc = \"SCS_DET_DSSC1M-1\"\n",
" dinstance = \"DSSC1M1\"\n",
"print(\"Detector in use is {}\".format(loc)) \n",
"\n",
"offset_image = int(offset_image)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-21T11:30:07.263445Z",
"start_time": "2019-02-21T11:30:07.217070Z"
},
"collapsed": true
},
"outputs": [],
"source": [
"def combine_stack(d, sdim):\n",
" combined = np.zeros((sdim, 2048,2048))\n",
" combined[...] = np.nan\n",
" d = np.moveaxis(d, 2, 3)\n",
" dy = 0\n",
" quad_pos = [\n",
" (0, 145),\n",
" (-5, 25),\n",
" (130, 15),\n",
" (130, 140),\n",
" ]\n",
" \n",
" px = 0.236\n",
" py = 0.204\n",
" with h5py.File(geo_file, \"r\") as gf:\n",
" \n",
" for i in range(16):\n",
" t1 = gf[\"Q{}M{}\"]\n",
"\n",
" if True: #try:\n",
" if i < 8:\n",
" dx = -512\n",
" if i > 3:\n",
" dx -= 25\n",
" mx = 1\n",
" my = i % 8\n",
" combined[:, my*128+dy:(my+1)*128+dy,\n",
" mx*512-dx:(mx+1)*512-dx] = np.rollaxis(d[i],2,1)[:,:,::-1]\n",
" dy += 30\n",
" if i == 3:\n",
" dy += 30\n",
"\n",
" elif i < 12:\n",
" dx = 80 - 50\n",
" if i == 8:\n",
" dy = 4*30 + 30 +50 -30\n",
"\n",
" mx = 1\n",
" my = i % 8 +4\n",
" combined[:, my*128+dy:(my+1)*128+dy,\n",
" mx*512-dx:(mx+1)*512-dx] = np.rollaxis(d[i],2,1)[:,::-1,:]\n",
" dy += 30\n",
" else:\n",
" dx = 100 - 50\n",
" if i == 11:\n",
" dy = 20\n",
"\n",
" mx = 1\n",
" my = i - 14\n",
"\n",
" combined[:, my*128+dy:(my+1)*128+dy,\n",
" mx*512-dx:(mx+1)*512-dx] = np.rollaxis(d[i],2,1)[:,::-1,:]\n",
" dy += 30\n",
" #except:\n",
"\n",
" # continue\n",
"\n",
" \n",
" return combined"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-21T11:30:07.974174Z",
"start_time": "2019-02-21T11:30:07.914832Z"
},
"collapsed": true
},
"outputs": [],
"source": [
"# set everything up filewise\n",
"from queue import Queue\n",
"from collections import OrderedDict\n",
"\n",
"def map_modules_from_files(filelist):\n",
" module_files = OrderedDict()\n",
" mod_ids = OrderedDict()\n",
" total_sequences = 0\n",
" sequences_qm = {}\n",
" one_module = None\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",
" sequences_qm[name] = 0\n",
" for file in filelist:\n",
" if file_infix in file:\n",
" if not one_module:\n",
" one_module = file, num\n",
" module_files[name].put(file)\n",
" total_sequences += 1\n",
" sequences_qm[name] += 1\n",
" \n",
" return module_files, mod_ids, total_sequences, sequences_qm, one_module\n",
"\n",
"dirlist = sorted(os.listdir(in_folder))\n",
"file_list = []\n",
"\n",
"\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, total_sequences, sequences_qm, one_module = map_modules_from_files(file_list)\n",
"MAX_PAR = min(MAX_PAR, total_sequences)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Processed Files ##"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-21T11:30:08.870802Z",
"start_time": "2019-02-21T11:30:08.826285Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"import copy\n",
"from IPython.display import HTML, display, Markdown, Latex\n",
"import tabulate\n",
"print(\"Processing a total of {} sequence files in chunks of {}\".format(total_sequences, MAX_PAR))\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",
"# restore the queue\n",
"mapped_files, mod_ids, total_sequences, sequences_qm, one_module = map_modules_from_files(file_list)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-21T11:30:16.057429Z",
"start_time": "2019-02-21T11:30:10.082114Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"import copy\n",
"from functools import partial\n",
"def correct_module(total_sequences, sequences_qm, loc, dinstance, offset_image,\n",
" mask_noisy_asic, mask_cold_asic, noisy_pix_threshold, chunksize, inp):\n",
" import numpy as np\n",
" import copy\n",
" import h5py\n",
" from cal_tools.enums import BadPixels\n",
" \n",
" filename, filename_out, channel, qm = inp\n",
" h5path = \"INSTRUMENT/{}/DET/{}CH0:xtdf/\".format(loc, channel)\n",
" h5path_idx = \"INDEX/{}/DET/{}CH0:xtdf/\".format(loc, channel)\n",
" \n",
" low_edges = None\n",
" hists_signal_low = None\n",
" high_edges = None\n",
" hists_signal_high = None\n",
" pulse_edges = None\n",
" \n",
" def copy_and_sanitize_non_cal_data(infile, outfile):\n",
" # these are touched in the correct function, do not copy them here\n",
" dont_copy = [\"data\"]\n",
" dont_copy = [h5path + \"image/{}\".format(do)\n",
" for do in dont_copy]\n",
"\n",
" # a visitor to copy everything else\n",
" def visitor(k, item):\n",
" if k not in dont_copy:\n",
"\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)\n",
"\n",
" try:\n",
" with h5py.File(filename, \"r\", driver=\"core\") as infile:\n",
" with h5py.File(filename_out, \"w\") as outfile:\n",
" copy_and_sanitize_non_cal_data(infile, outfile)\n",
" # get indices of last images in each train\n",
" first_arr = np.squeeze(infile[h5path_idx+\"image/first\"]).astype(np.int)\n",
" last_arr = np.concatenate((first_arr[1:], np.array([-1,]))).astype(np.int)\n",
" assert first_arr.size == last_arr.size\n",
" oshape = list(infile[h5path+\"image/data\"].shape)\n",
" if len(oshape) == 4:\n",
" oshape = [oshape[0],]+oshape[2:]\n",
" chunks = (chunksize, oshape[1], oshape[2])\n",
" ddset = outfile.create_dataset(h5path + \"image/data\",\n",
" oshape, chunks=chunks,\n",
" dtype=np.uint16,\n",
" fletcher32=True)\n",
" \n",
" mdset = outfile.create_dataset(h5path + \"image/mask\",\n",
" oshape, chunks=chunks,\n",
" dtype=np.uint32,\n",
" compression=\"gzip\",\n",
" compression_opts=1,\n",
" shuffle=True,\n",
" fletcher32=True)\n",
" \n",
" for train in range(first_arr.size):\n",
" first = first_arr[train]\n",
" last = last_arr[train]\n",
" data = np.squeeze(infile[h5path+\"image/data\"][first:last, ...].astype(np.float32))\n",
" pulseId = np.squeeze(infile[h5path+\"image/pulseId\"][first:last, ...])\n",
" data -= data[offset_image, ...]\n",
" \n",
" if train == 0:\n",
" pulseId = np.repeat(pulseId[:, None], data.shape[1], axis=1)\n",
" pulseId = np.repeat(pulseId[:,:,None], data.shape[2], axis=2)\n",
" bins = (55, pulseId.max())\n",
" rnge = [[-5, 50], [0, pulseId.max()]]\n",
" \n",
" hists_signal_low, low_edges, pulse_edges = np.histogram2d(data.flatten(),\n",
" pulseId.flatten(),\n",
" bins=bins,\n",
" range=rnge)\n",
" rnge = [[-5, 300], [0, pulseId.max()]]\n",
" hists_signal_high, high_edges, _ = np.histogram2d(data.flatten(),\n",
" pulseId.flatten(),\n",
" bins=bins,\n",
" range=rnge) \n",
" data[data < 0] = 0\n",
" ddset[first:last, ...] = data.astype(np.uint16)\n",
" \n",
" # find static and noisy values in dark images\n",
" data = infile[h5path+\"image/data\"][last, ...].astype(np.float32)\n",
" bpix = np.zeros(oshape[1:], np.uint32)\n",
" dark_std = np.std(data, axis=0)\n",
" bpix[dark_std > noisy_pix_threshold] = BadPixels.NOISE_OUT_OF_THRESHOLD.value\n",
" \n",
" for i in range(8):\n",
" for j in range(2):\n",
" count_noise = np.count_nonzero(bpix[i*64:(i+1)*64, j*64:(j+1)*64])\n",
" asic_std = np.std(data[:, i*64:(i+1)*64, j*64:(j+1)*64])\n",
" if mask_noisy_asic:\n",
" if count_noise/(64*64) > mask_noisy_asic:\n",
" bpix[i*64:(i+1)*64, j*64:(j+1)*64] = BadPixels.NOISY_ADC.value\n",
"\n",
" if mask_cold_asic:\n",
" count_cold = np.count_nonzero(asic_std < 0.5)\n",
" if count_cold/(64*64) > mask_cold_asic:\n",
" bpix[i*64:(i+1)*64, j*64:(j+1)*64] = BadPixels.ASIC_STD_BELOW_NOISE.value\n",
" \n",
" mdset[...] = np.repeat(bpix[None,...], infile[h5path+\"image/data\"].shape[0], axis=0)\n",
" \n",
" except Exception as e:\n",
" print(e)\n",
" success = False\n",
" reason = \"Error\"\n",
" \n",
" \n",
" \n",
" return (hists_signal_low, hists_signal_high, low_edges, high_edges, pulse_edges)\n",
" \n",
"done = False\n",
"first_files = []\n",
"inp = []\n",
"left = total_sequences\n",
"\n",
"hists_signal_low = 0\n",
"hists_signal_high = 0 \n",
"\n",
"low_edges, high_edges, pulse_edges = None, None, None\n",
"\n",
"\n",
"while not done:\n",
" \n",
" dones = []\n",
" first = True\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",
" dones.append(mapped_files[qm].empty())\n",
" else:\n",
" print(\"Skipping {}\".format(qm))\n",
" first_files.append((None, None))\n",
" continue\n",
" fout = os.path.abspath(\"{}/{}\".format(out_folder, (os.path.split(fname_in)[-1]).replace(\"RAW\", \"CORR\")))\n",
" if first:\n",
" first_files.append((fname_in, fout))\n",
" inp.append((fname_in, fout, i, qm))\n",
" first = False\n",
" if len(inp) >= min(MAX_PAR, left):\n",
" print(\"Running {} tasks parallel\".format(len(inp)))\n",
" p = partial(correct_module, total_sequences, sequences_qm,\n",
" loc, dinstance, offset_image, mask_noisy_asic,\n",
" mask_cold_asic, noisy_pix_threshold, chunk_size_idim)\n",
" \n",
" r = view.map_sync(p, inp)\n",
" #r = list(map(p, inp))\n",
" inp = []\n",
" left -= MAX_PAR\n",
" \n",
" for rr in r:\n",
" if rr is not None:\n",
" hl, hh, low_edges, high_edges, pulse_edges = rr \n",
" if hl is not None: # any one being None will also make the others None\n",
" hists_signal_low += hl.astype(np.float64)\n",
" hists_signal_high += hh.astype(np.float64) \n",
" \n",
" done = all(dones)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:28:51.765030Z",
"start_time": "2019-02-18T17:28:51.714783Z"
},
"collapsed": true
},
"outputs": [],
"source": [
"from mpl_toolkits.mplot3d import Axes3D\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib import cm\n",
"from matplotlib.ticker import LinearLocator, FormatStrFormatter\n",
"import numpy as np\n",
"%matplotlib inline\n",
"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",
"\n",
" Z = data.T\n",
"\n",
" # Plot the surface.\n",
" surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,\n",
" linewidth=0, antialiased=False)\n",
" ax.set_xlabel(x_axis)\n",
" ax.set_ylabel(y_axis)\n",
" ax.set_zlabel(\"Counts\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:28:53.690522Z",
"start_time": "2019-02-18T17:28:52.860143Z"
},
"collapsed": true
},
"outputs": [],
"source": [
"def do_2d_plot(data, edges, y_axis, x_axis):\n",
" from matplotlib.colors import LogNorm\n",
" fig = plt.figure(figsize=(10,10))\n",
" ax = fig.add_subplot(111)\n",
" extent = [np.min(edges[1]), np.max(edges[1]),np.min(edges[0]), np.max(edges[0])]\n",
" im = ax.imshow(data[::-1,:], extent=extent, aspect=\"auto\", norm=LogNorm(vmin=1, vmax=np.max(data)))\n",
" ax.set_xlabel(x_axis)\n",
" ax.set_ylabel(y_axis)\n",
" cb = fig.colorbar(im)\n",
" cb.set_label(\"Counts\")\n",
" \n",
" \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mean Intensity per Pulse ##\n",
"\n",
"The following plots show the mean signal for each pulse in a detailed and expanded intensity region."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:28:57.327702Z",
"start_time": "2019-02-18T17:28:54.377061Z"
},
"collapsed": false,
"scrolled": false
},
"outputs": [],
"source": [
"do_3d_plot(hists_signal_low, [low_edges, pulse_edges], \"Signal (ADU)\", \"Pulse id\")\n",
"do_2d_plot(hists_signal_low, [low_edges, pulse_edges], \"Signal (ADU)\", \"Pulse id\")\n",
"do_3d_plot(hists_signal_high, [high_edges, pulse_edges], \"Signal (ADU)\", \"Pulse id\")\n",
"do_2d_plot(hists_signal_high, [high_edges, pulse_edges], \"Signal (ADU)\", \"Pulse id\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:29:20.634480Z",
"start_time": "2019-02-18T17:28:57.329231Z"
},
"collapsed": true
},
"outputs": [],
"source": [
"\n",
"corrected = []\n",
"raw = []\n",
"mask = []\n",
"cell_fac = 1\n",
"first_idx = 400*10+40\n",
"last_idx = 400*10+56\n",
"pulse_ids = []\n",
"train_ids = []\n",
"for i, ff in enumerate(first_files[:16]):\n",
" try:\n",
"\n",
" rf, cf = ff\n",
" #print(cf, i)\n",
" if rf is None:\n",
" \n",
" raise Exception(\"File not present\")\n",
" infile = h5py.File(rf, \"r\")\n",
" raw.append(np.array(infile[\"/INSTRUMENT/{}/DET/{}CH0:xtdf/image/data\".format(loc, i)][first_idx:last_idx,0,...]))\n",
" infile.close()\n",
" \n",
" infile = h5py.File(cf, \"r\")\n",
" corrected.append(np.array(infile[\"/INSTRUMENT/{}/DET/{}CH0:xtdf/image/data\".format(loc, i)][first_idx:last_idx,...]))\n",
" mask.append(np.array(infile[\"/INSTRUMENT/{}/DET/{}CH0:xtdf/image/mask\".format(loc, i)][first_idx:last_idx,...]))\n",
" pulse_ids.append(np.squeeze(infile[\"/INSTRUMENT/{}/DET/{}CH0:xtdf/image/pulseId\".format(loc, i)][first_idx:last_idx,...]))\n",
" train_ids.append(np.squeeze(infile[\"/INSTRUMENT/{}/DET/{}CH0:xtdf/image/trainId\".format(loc, i)][first_idx:last_idx,...]))\n",
" infile.close()\n",
" \n",
" except Exception as e:\n",
" print(e)\n",
" corrected.append(np.zeros((last_idx-first_idx, 512, 128)))\n",
" mask.append(np.zeros((last_idx-first_idx, 512, 128)))\n",
" raw.append(np.zeros((last_idx-first_idx, 512, 128)))\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def combine_stack(d, sdim):\n",
" combined = np.zeros((sdim, 1300,1300), np.float32)\n",
" combined[...] = 0#np.nan\n",
" \n",
" dy = 0\n",
" quad_pos = [\n",
" (0, 145),\n",
" (130, 140),\n",
" (125, 15),\n",
" (0, 15),\n",
" \n",
" ]\n",
" \n",
" px = 0.236\n",
" py = 0.204\n",
" with h5py.File(geo_file, \"r\") as gf:\n",
" \n",
" for i in range(16):\n",
" mi = i\n",
" #if i // 4 == 0 or i // 4 == 1:\n",
" mi = 3-(i%4)\n",
" mp = gf[\"Q{}/M{}/Position\".format(i//4+1, mi%4+1)][()]\n",
" t1 = gf[\"Q{}/M{}/T01/Position\".format(i//4+1, i%4+1)][()]\n",
" t2 = gf[\"Q{}/M{}/T02/Position\".format(i//4+1, i%4+1)][()]\n",
" if i//4 < 2:\n",
" t1, t2 = t2, t1\n",
" \n",
" if i // 4 == 0 or i // 4 == 1:\n",
" td = d[i][:,::-1,:]\n",
" else:\n",
" td = d[i][:,:,::-1]\n",
" \n",
" t1d = td[:,:,:256]\n",
" t2d = td[:,:,256:]\n",
" \n",
" x0t1 = int((t1[0]+mp[0])/px)\n",
" y0t1 = int((t1[1]+mp[1])/py)\n",
" x0t2 = int((t2[0]+mp[0])/px)\n",
" y0t2 = int((t2[1]+mp[1])/py)\n",
" \n",
" x0t1 += int(quad_pos[i//4][1]/px)\n",
" x0t2 += int(quad_pos[i//4][1]/px)\n",
" y0t1 += int(quad_pos[i//4][0]/py)+combined.shape[1]//16\n",
" y0t2 += int(quad_pos[i//4][0]/py)+combined.shape[1]//16\n",
" combined[:,y0t1:y0t1+128,x0t1:x0t1+256] = t1d\n",
" combined[:,y0t2:y0t2+128,x0t2:x0t2+256] = t2d\n",
"\n",
" return combined\n",
"\n",
"combined = combine_stack(corrected, last_idx-first_idx)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:29:27.025667Z",
"start_time": "2019-02-18T17:29:20.642029Z"
},
"collapsed": true
},
"outputs": [],
"source": [
"\n",
"combined = combine_stack(corrected, last_idx-first_idx)\n",
"combined_raw = combine_stack(raw, last_idx-first_idx)\n",
"combined_mask = combine_stack(mask, last_idx-first_idx)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mean RAW Preview ###\n",
"\n",
"The per pixel mean of the first 128 images of the RAW data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:29:33.226396Z",
"start_time": "2019-02-18T17:29:27.027758Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"fig = plt.figure(figsize=(20,10))\n",
"ax = fig.add_subplot(111)\n",
"im = ax.imshow(np.mean(combined_raw[:,...],axis=0),\n",
" vmin=min(0.75*np.median(combined_raw[combined_raw > 0]), -5),\n",
" vmax=max(1.5*np.median(combined_raw[combined_raw > 0]), 50), cmap=\"jet\")\n",
"cb = fig.colorbar(im, ax=ax)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Single Shot Preview ###\n",
"\n",
"A single shot image from cell 2 of the first train"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:29:33.761015Z",
"start_time": "2019-02-18T17:29:33.227922Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"\n",
"fig = plt.figure(figsize=(20,10))\n",
"ax = fig.add_subplot(111)\n",
"dim = combined[2,...]\n",
"\n",
"im = ax.imshow(dim, vmin=-0, vmax=max(1.5*np.median(dim[dim > 0]), 50), cmap=\"jet\", interpolation=\"nearest\")\n",
"cb = fig.colorbar(im, ax=ax)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:29:35.903487Z",
"start_time": "2019-02-18T17:29:33.762568Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20,10))\n",
"ax = fig.add_subplot(111)\n",
"h = ax.hist(dim.flatten(), bins=100, range=(0, 100))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mean CORRECTED Preview ###\n",
"\n",
"The per pixel mean of the first 128 images of the CORRECTED data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:29:39.369686Z",
"start_time": "2019-02-18T17:29:35.905152Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"\n",
"fig = plt.figure(figsize=(20,10))\n",
"ax = fig.add_subplot(111)\n",
"im = ax.imshow(np.mean(combined[:,...], axis=0), vmin=0,\n",
" vmax=max(1.5*np.median(combined[combined > 0]), 10), cmap=\"jet\", interpolation=\"nearest\")\n",
"cb = fig.colorbar(im, ax=ax)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Max CORRECTED Preview ###\n",
"\n",
"The per pixel maximum of the first 128 images of the CORRECTED data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"\n",
"fig = plt.figure(figsize=(20,10))\n",
"ax = fig.add_subplot(111)\n",
"im = ax.imshow(np.max(combined[:,...], axis=0), vmin=0,\n",
" vmax=max(100*np.median(combined[combined > 0]), 20), cmap=\"jet\", interpolation=\"nearest\")\n",
"cb = fig.colorbar(im, ax=ax)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:29:49.217848Z",
"start_time": "2019-02-18T17:29:39.371232Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20,10))\n",
"ax = fig.add_subplot(111)\n",
"combined[combined <= 0] = 0\n",
"h = ax.hist(combined.flatten(), bins=100, range=(-5, 100), log=True)\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Bad Pixels ##\n",
"The mask contains dedicated entries for all pixels and memory cells as well as all three gains stages. Each mask entry is encoded in 32 bits as:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:29:49.651913Z",
"start_time": "2019-02-18T17:29:49.643556Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"from cal_tools.enums import BadPixels\n",
"from IPython.display import HTML, display, Markdown, Latex\n",
"import tabulate\n",
"table = []\n",
"for item in BadPixels:\n",
" table.append((item.name, \"{:016b}\".format(item.value)))\n",
"md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=[\"Bad pixel type\", \"Bit mask\"])))"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### Full Train Bad Pixels ###"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:29:51.686562Z",
"start_time": "2019-02-18T17:29:50.088883Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20,10))\n",
"ax = fig.add_subplot(111)\n",
"im = ax.imshow(np.log2(np.max(combined_mask[:,...], axis=0)), vmin=0,\n",
" vmax=32, cmap=\"jet\")\n",
"cb = fig.colorbar(im, ax=ax)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Full Train Bad Pixels - Only Dark Char. Related ###"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-02-18T17:29:53.662423Z",
"start_time": "2019-02-18T17:29:51.688376Z"
},
"collapsed": false
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(20,10))\n",
"ax = fig.add_subplot(111)\n",
"im = ax.imshow(np.max((combined_mask.astype(np.uint32)[:,...] & BadPixels.NOISY_ADC.value) != 0, axis=0), vmin=0,\n",
" vmax=1, cmap=\"jet\")\n",
"cb = fig.colorbar(im, ax=ax)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"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": 2
}