From a59747f3f0f41d65466a8076c6fe09327343526a Mon Sep 17 00:00:00 2001
From: Steffen Hauf <steffen.hauf@xfel.eu>
Date: Thu, 27 Jun 2019 16:36:44 +0200
Subject: [PATCH] Propagate backlog of DSSC notebook related changes from
production system as of 06/2019
---
.../DSSC/Characterize_DSSC_Darks_NBC.ipynb | 573 ++++++++++
notebooks/DSSC/DSSC_Correct_and_Verify.ipynb | 1004 +++++++++++++++++
2 files changed, 1577 insertions(+)
create mode 100644 notebooks/DSSC/Characterize_DSSC_Darks_NBC.ipynb
create mode 100644 notebooks/DSSC/DSSC_Correct_and_Verify.ipynb
diff --git a/notebooks/DSSC/Characterize_DSSC_Darks_NBC.ipynb b/notebooks/DSSC/Characterize_DSSC_Darks_NBC.ipynb
new file mode 100644
index 000000000..f6d083c2f
--- /dev/null
+++ b/notebooks/DSSC/Characterize_DSSC_Darks_NBC.ipynb
@@ -0,0 +1,573 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Characterize Dark Images #\n",
+ "\n",
+ "Author: S. Hauf, Version: 0.1\n",
+ "\n",
+ "The following code analyzes a set of dark images taken with the AGIPD detector to deduce detector offsets and noise. Data for the detector's three gain stages needs to be present, separated into separate runs.\n",
+ "\n",
+ "The notebook explicitely does what pyDetLib provides in its offset calculation method for streaming data."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-02-20T12:42:51.255184Z",
+ "start_time": "2019-02-20T12:42:51.225500Z"
+ },
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "cluster_profile = \"noDB\" # The ipcluster profile to use\n",
+ "in_folder = \"/gpfs/exfel/exp/SCS/201930/p900079/raw\" # path to input data, required\n",
+ "out_folder = \"/gpfs/exfel/data/scratch/haufs/test/\" # path to output to, required\n",
+ "sequences = [0] # sequence files to evaluate.\n",
+ "\n",
+ "run = 20 # run number in which data was recorded, required\n",
+ "\n",
+ "mem_cells = 0 # number of memory cells used, set to 0 to automatically infer\n",
+ "local_output = False # output constants locally\n",
+ "db_output = True # output constants to database\n",
+ "bias_voltage = 300 # detector bias voltage\n",
+ "cal_db_interface = \"tcp://max-exfl016:8020\" # the database interface to use\n",
+ "rawversion = 2 # RAW file format version\n",
+ "dont_use_dir_date = False # don't use the dir creation date for determining the creation time\n",
+ "\n",
+ "thresholds_offset_sigma = 3. # thresholds in terms of n sigma noise for offset deduced bad pixels\n",
+ "thresholds_offset_hard = [4000, 8500] # thresholds in absolute ADU terms for offset deduced bad pixels\n",
+ "\n",
+ "thresholds_noise_sigma = 5. # thresholds in terms of n sigma noise for offset deduced bad pixels\n",
+ "thresholds_noise_hard = [4, 20] # thresholds in absolute ADU terms for offset deduced bad pixels\n",
+ "\n",
+ "instrument = \"SCS\" # the instrument\n",
+ "high_res_badpix_3d = False # set this to True if you need high-resolution 3d bad pixel plots. Runtime: ~ 1h\n",
+ "modules = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] # module to run for"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-02-20T12:42:52.599660Z",
+ "start_time": "2019-02-20T12:42:51.472138Z"
+ },
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# imports and things that do not usually need to be changed\n",
+ "from datetime import datetime\n",
+ "import warnings\n",
+ "warnings.filterwarnings('ignore')\n",
+ "from collections import OrderedDict\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",
+ "%matplotlib inline\n",
+ "\n",
+ "from cal_tools.tools import gain_map_files, parse_runs, run_prop_seq_from_path, get_notebook_name, get_dir_creation_date\n",
+ "from cal_tools.influx import InfluxLogger\n",
+ "from cal_tools.enums import BadPixels\n",
+ "from cal_tools.plotting import show_overview, plot_badpix_3d, create_constant_overview\n",
+ "\n",
+ "# make sure a cluster is running with ipcluster start --n=32, give it a while to start\n",
+ "from ipyparallel import Client\n",
+ "\n",
+ "view = Client(profile=cluster_profile)[:]\n",
+ "view.use_dill()\n",
+ "\n",
+ "from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions\n",
+ "\n",
+ "\n",
+ "# no need to change this\n",
+ "\n",
+ "QUADRANTS = 4\n",
+ "MODULES_PER_QUAD = 4\n",
+ "DET_FILE_INSET = \"DSSC\"\n",
+ "\n",
+ "max_cells = mem_cells\n",
+ " \n",
+ "offset_runs = OrderedDict()\n",
+ "offset_runs[\"high\"] = parse_runs(run)[0]\n",
+ "\n",
+ "creation_time=None\n",
+ "if not dont_use_dir_date:\n",
+ " creation_time = get_dir_creation_date(in_folder, run)\n",
+ "\n",
+ "\n",
+ "run, prop, seq = run_prop_seq_from_path(in_folder)\n",
+ "logger = InfluxLogger(detector=\"DSSC\", instrument=instrument, mem_cells=mem_cells,\n",
+ " notebook=get_notebook_name(), proposal=prop)\n",
+ "\n",
+ "print(\"Using {} as creation time of constant.\".format(creation_time))\n",
+ "\n",
+ "loc = None\n",
+ "if instrument == \"SCS\":\n",
+ " loc = \"SCS_DET_DSSC1M-1\"\n",
+ " dinstance = \"DSSC1M1\"\n",
+ "\n",
+ "print(\"Detector in use is {}\".format(loc)) "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-02-20T12:42:52.608214Z",
+ "start_time": "2019-02-20T12:42:52.601257Z"
+ },
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "print(\"Parameters are:\")\n",
+ "print(\"Proposal: {}\".format(prop))\n",
+ "print(\"Memory cells: {}/{}\".format(mem_cells, max_cells))\n",
+ "print(\"Runs: {}\".format([ v for v in offset_runs.values()]))\n",
+ "print(\"Sequences: {}\".format(sequences))\n",
+ "print(\"Using DB: {}\".format(db_output))\n",
+ "print(\"Input: {}\".format(in_folder))\n",
+ "print(\"Output: {}\".format(out_folder))\n",
+ "print(\"Bias voltage: {}V\".format(bias_voltage))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The following lines will create a queue of files which will the be executed module-parallel. Distiguishing between different gains."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-02-20T12:42:54.024731Z",
+ "start_time": "2019-02-20T12:42:53.901555Z"
+ },
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# set everything up filewise\n",
+ "if not os.path.exists(out_folder):\n",
+ " os.makedirs(out_folder)\n",
+ "\n",
+ "gmf = gain_map_files(in_folder, offset_runs, sequences, DET_FILE_INSET, QUADRANTS, MODULES_PER_QUAD)\n",
+ "gain_mapped_files, total_sequences, total_file_size = gmf\n",
+ "\n",
+ "print(\"Will process at total of {} sequences: {:0.2f} GB of data.\".format(total_sequences, total_file_size))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Calculate Offsets, Noise and Thresholds ##\n",
+ "\n",
+ "The calculation is performed per-pixel and per-memory-cell. Offsets are simply the median value for a set of dark data taken at a given gain, noise the standard deviation, and gain-bit values the medians of the gain array."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-02-20T10:50:55.839958Z",
+ "start_time": "2019-02-20T10:50:55.468134Z"
+ },
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "import copy\n",
+ "from functools import partial\n",
+ "def characterize_module(cells, bp_thresh, rawversion, loc, inp):\n",
+ " import numpy as np\n",
+ " import copy\n",
+ " import h5py\n",
+ " from cal_tools.enums import BadPixels\n",
+ " \n",
+ " def get_num_cells(fname, loc, module):\n",
+ " with h5py.File(fname, \"r\") as f:\n",
+ "\n",
+ " cells = f[\"INSTRUMENT/{}/DET/{}CH0:xtdf/image/cellId\".format(loc, module)][()]\n",
+ " maxcell = np.max(cells)\n",
+ " options = [100, 200, 400, 500, 600, 700, 800]\n",
+ " dists = [abs(o-maxcell) for o in options]\n",
+ " return options[np.argmin(dists)]\n",
+ " \n",
+ " filename, filename_out, channel = inp\n",
+ " \n",
+ "\n",
+ " if cells == 0:\n",
+ " cells = get_num_cells(filename, loc, channel)\n",
+ "\n",
+ " #print(\"Using {} memory cells\".format(cells))\n",
+ " \n",
+ " thresholds_offset_hard, thresholds_offset_sigma, thresholds_noise_hard, thresholds_noise_sigma = bp_thresh \n",
+ "\n",
+ " infile = h5py.File(filename, \"r\", driver=\"core\")\n",
+ " if rawversion == 2:\n",
+ " count = np.squeeze(infile[\"/INDEX/{}/DET/{}CH0:xtdf/image/count\".format(loc, channel)])\n",
+ " first = np.squeeze(infile[\"/INDEX/{}/DET/{}CH0:xtdf/image/first\".format(loc, channel)])\n",
+ " last_index = int(first[count != 0][-1]+count[count != 0][-1])\n",
+ " first_index = int(first[count != 0][0])\n",
+ " else:\n",
+ " status = np.squeeze(infile[\"/INDEX/{}/DET/{}CH0:xtdf/image/status\".format(loc, channel)])\n",
+ " if np.count_nonzero(status != 0) == 0:\n",
+ " return\n",
+ " last = np.squeeze(infile[\"/INDEX/{}/DET/{}CH0:xtdf/image/last\".format(loc, channel)])\n",
+ " first = np.squeeze(infile[\"/INDEX/{}/DET/{}CH0:xtdf/image/first\".format(loc, channel)])\n",
+ " last_index = int(last[status != 0][-1]) + 1\n",
+ " first_index = int(first[status != 0][0])\n",
+ " im = np.array(infile[\"/INSTRUMENT/{}/DET/{}CH0:xtdf/image/data\".format(loc, channel)][first_index:last_index,...]) \n",
+ " cellIds = np.squeeze(infile[\"/INSTRUMENT/{}/DET/{}CH0:xtdf/image/cellId\".format(loc, channel)][first_index:last_index,...]) \n",
+ " \n",
+ " infile.close()\n",
+ "\n",
+ " \n",
+ " im = im[:, 0, ...].astype(np.float32)\n",
+ " \n",
+ " im = np.rollaxis(im, 2)\n",
+ " im = np.rollaxis(im, 2, 1)\n",
+ "\n",
+ " mcells = cells #max(cells, np.max(cellIds)+1)\n",
+ " offset = np.zeros((im.shape[0], im.shape[1], mcells))\n",
+ " noise = np.zeros((im.shape[0], im.shape[1], mcells))\n",
+ " \n",
+ " for cc in np.unique(cellIds[cellIds < mcells]):\n",
+ " cellidx = cellIds == cc\n",
+ " offset[...,cc] = np.median(im[..., cellidx], axis=2)\n",
+ " noise[...,cc] = np.std(im[..., cellidx], axis=2)\n",
+ " \n",
+ " \n",
+ " # bad pixels\n",
+ " bp = np.zeros(offset.shape, np.uint32)\n",
+ " # offset related bad pixels\n",
+ " offset_mn = np.nanmedian(offset, axis=(0,1))\n",
+ " offset_std = np.nanstd(offset, axis=(0,1)) \n",
+ " \n",
+ " bp[(offset < offset_mn-thresholds_offset_sigma*offset_std) |\n",
+ " (offset > offset_mn+thresholds_offset_sigma*offset_std)] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value\n",
+ " bp[(offset < thresholds_offset_hard[0]) | (offset > thresholds_offset_hard[1])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value\n",
+ " bp[~np.isfinite(offset)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value\n",
+ " \n",
+ " # noise related bad pixels\n",
+ " noise_mn = np.nanmedian(noise, axis=(0,1))\n",
+ " noise_std = np.nanstd(noise, axis=(0,1)) \n",
+ " \n",
+ " bp[(noise < noise_mn-thresholds_noise_sigma*noise_std) |\n",
+ " (noise > noise_mn+thresholds_noise_sigma*noise_std)] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value\n",
+ " bp[(noise < thresholds_noise_hard[0]) | (noise > thresholds_noise_hard[1])] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value\n",
+ " bp[~np.isfinite(noise)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value\n",
+ "\n",
+ "\n",
+ " return offset, noise, bp, cells\n",
+ " \n",
+ " \n",
+ "offset_g = OrderedDict()\n",
+ "noise_g = OrderedDict()\n",
+ "gain_g = OrderedDict()\n",
+ "badpix_g = OrderedDict()\n",
+ "gg = 0\n",
+ "\n",
+ "start = datetime.now()\n",
+ "all_cells = []\n",
+ "\n",
+ "for gain, mapped_files in gain_mapped_files.items():\n",
+ " \n",
+ " inp = []\n",
+ " dones = []\n",
+ " for i in modules:\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 = mapped_files[qm].get() \n",
+ " dones.append(mapped_files[qm].empty())\n",
+ " else:\n",
+ " continue\n",
+ " fout = os.path.abspath(\"{}/{}\".format(out_folder, (os.path.split(fname_in)[-1]).replace(\"RAW\", \"CORR\")))\n",
+ " inp.append((fname_in, fout, i))\n",
+ " first = False\n",
+ " p = partial(characterize_module, max_cells,\n",
+ " (thresholds_offset_hard, thresholds_offset_sigma,\n",
+ " thresholds_noise_hard, thresholds_noise_sigma), rawversion, loc)\n",
+ " results = list(map(p, inp))\n",
+ " #results = view.map_sync(p, inp)\n",
+ " for ii, r in enumerate(results):\n",
+ " i = modules[ii]\n",
+ " offset, noise, bp, thiscell = r\n",
+ " all_cells.append(thiscell)\n",
+ " qm = \"Q{}M{}\".format(i//4 +1, i % 4 + 1)\n",
+ " if qm not in offset_g:\n",
+ " offset_g[qm] = np.zeros((offset.shape[0], offset.shape[1], offset.shape[2]))\n",
+ " noise_g[qm] = np.zeros_like(offset_g[qm])\n",
+ " \n",
+ " badpix_g[qm] = np.zeros_like(offset_g[qm], np.uint32)\n",
+ " \n",
+ " offset_g[qm][...] = offset\n",
+ " noise_g[qm][...] = noise\n",
+ " badpix_g[qm][...] = bp\n",
+ " gg +=1\n",
+ "\n",
+ "duration = (datetime.now()-start).total_seconds()\n",
+ "logger.runtime_summary_entry(success=True, runtime=duration,\n",
+ " total_sequences=total_sequences,\n",
+ " filesize=total_file_size)\n",
+ "logger.send()\n",
+ "max_cells = np.max(all_cells)\n",
+ "print(\"Using {} memory cells\".format(max_cells))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2018-12-06T09:38:18.234582Z",
+ "start_time": "2018-12-06T09:38:18.222838Z"
+ },
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "res = OrderedDict()\n",
+ "for i in modules:\n",
+ " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
+ " res[qm] = {'Offset': offset_g[qm],\n",
+ " 'Noise': noise_g[qm],\n",
+ " 'BadPixels': badpix_g[qm] \n",
+ " }\n",
+ " \n",
+ "if local_output:\n",
+ " for qm in offset_g.keys():\n",
+ " ofile = \"{}/dssc_offset_store_{}_{}.h5\".format(out_folder, \"_\".join(offset_runs.values()), qm)\n",
+ " store_file = h5py.File(ofile, \"w\")\n",
+ " store_file[\"{}/Offset/0/data\".format(qm)] = offset_g[qm]\n",
+ " store_file[\"{}/Noise/0/data\".format(qm)] = noise_g[qm]\n",
+ " store_file[\"{}/BadPixels/0/data\".format(qm)] = badpix_g[qm]\n",
+ " store_file.close()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2018-12-06T09:49:32.449330Z",
+ "start_time": "2018-12-06T09:49:20.231607Z"
+ },
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "if db_output:\n",
+ " for qm in offset_g.keys():\n",
+ " metadata = ConstantMetaData()\n",
+ " offset = Constants.DSSC.Offset()\n",
+ " offset.data = offset_g[qm]\n",
+ " metadata.calibration_constant = offset\n",
+ "\n",
+ " # set the operating condition\n",
+ " condition = Conditions.Dark.DSSC(memory_cells=max_cells, bias_voltage=bias_voltage)\n",
+ " detinst = getattr(Detectors, dinstance)\n",
+ "\n",
+ " device = getattr(detinst, qm)\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, start=creation_time)\n",
+ " metadata.send(cal_db_interface, timeout=3000000)\n",
+ " \n",
+ " \n",
+ " metadata = ConstantMetaData()\n",
+ " noise = Constants.DSSC.Noise()\n",
+ " noise.data = noise_g[qm]\n",
+ " metadata.calibration_constant = noise\n",
+ "\n",
+ " # set the operating condition\n",
+ " condition = Conditions.Dark.DSSC(memory_cells=max_cells, bias_voltage=bias_voltage)\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, start=creation_time)\n",
+ " metadata.send(cal_db_interface, timeout=3000000)\n",
+ " \n",
+ " continue # no bad pixels yet\n",
+ " metadata = ConstantMetaData()\n",
+ " badpixels = Constants.DSSC.BadPixelsDark()\n",
+ " badpixels.data = badpix_g[qm]\n",
+ " metadata.calibration_constant = badpixels\n",
+ "\n",
+ " # set the operating condition\n",
+ " condition = Conditions.Dark.DSSC(memory_cells=max_cells, bias_voltage=bias_voltage)\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, start=creation_time)\n",
+ " metadata.send(cal_db_interface, timeout=3000000)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Single-Cell Overviews ##\n",
+ "\n",
+ "Single cell overviews allow to identify potential effects on all memory cells, e.g. on sensor level. Additionally, they should serve as a first sanity check on expected behaviour, e.g. if structuring on the ASIC level is visible in the offsets, but otherwise no immediate artifacts are visible."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2018-12-06T09:49:14.540552Z",
+ "start_time": "2018-12-06T09:49:13.009674Z"
+ },
+ "collapsed": false,
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "cell = 3\n",
+ "gain = 0\n",
+ "out_folder = None\n",
+ "show_overview(res, cell, gain, out_folder=out_folder, infix=\"_\".join(offset_runs.values()))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Global Bad Pixel Behaviour ##\n",
+ "\n",
+ "The following plots show the results of bad pixel evaluation for all evaluated memory cells. Cells are stacked in the Z-dimension, while pixels values in x/y are rebinned with a factor of 2. This excludes single bad pixels present only in disconnected pixels. Hence, any bad pixels spanning at least 4 pixels in the x/y-plane, or across at least two memory cells are indicated. Colors encode the bad pixel type, or mixed type."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true,
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "cols = {BadPixels.NOISE_OUT_OF_THRESHOLD.value: (BadPixels.NOISE_OUT_OF_THRESHOLD.name, '#FF000080'),\n",
+ " BadPixels.OFFSET_NOISE_EVAL_ERROR.value: (BadPixels.OFFSET_NOISE_EVAL_ERROR.name, '#0000FF80'),\n",
+ " BadPixels.OFFSET_OUT_OF_THRESHOLD.value: (BadPixels.OFFSET_OUT_OF_THRESHOLD.name, '#00FF0080'),\n",
+ " BadPixels.OFFSET_OUT_OF_THRESHOLD.value | BadPixels.NOISE_OUT_OF_THRESHOLD.value: ('MIXED', '#DD00DD80')}\n",
+ "\n",
+ "rebin = 8 if not high_res_badpix_3d else 2\n",
+ "\n",
+ "gain = 0\n",
+ "for mod, data in badpix_g.items():\n",
+ " plot_badpix_3d(data[...,gain], cols, title=mod, rebin_fac=rebin)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Aggregate values, and per Cell behaviour ##\n",
+ "\n",
+ "The following tables and plots give an overview of statistical aggregates for each constant, as well as per cell behavior."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true,
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "create_constant_overview(offset_g, \"Offset (ADU)\", max_cells, 4000, 8000,\n",
+ " out_folder=out_folder, infix=\"_\".join(offset_runs.values()))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true,
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "create_constant_overview(noise_g, \"Noise (ADU)\", max_cells, 0, 100,\n",
+ " out_folder=out_folder, infix=\"_\".join(offset_runs.values()))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "bad_pixel_aggregate_g = OrderedDict()\n",
+ "for m, d in badpix_g.items():\n",
+ " bad_pixel_aggregate_g[m] = d.astype(np.bool).astype(np.float)\n",
+ "create_constant_overview(bad_pixel_aggregate_g, \"Bad pixel fraction\", max_cells, 0, 0.10, 3,\n",
+ " out_folder=out_folder, infix=\"_\".join(offset_runs.values()))"
+ ]
+ },
+ {
+ "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
+}
diff --git a/notebooks/DSSC/DSSC_Correct_and_Verify.ipynb b/notebooks/DSSC/DSSC_Correct_and_Verify.ipynb
new file mode 100644
index 000000000..afa707fb4
--- /dev/null
+++ b/notebooks/DSSC/DSSC_Correct_and_Verify.ipynb
@@ -0,0 +1,1004 @@
+{
+ "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
+}
--
GitLab