From 6b137177865be46f2d6addf3a203c345da4a4ad7 Mon Sep 17 00:00:00 2001 From: Karim Ahmed <karim.ahmed@xfel.eu> Date: Fri, 27 Sep 2019 15:58:17 +0200 Subject: [PATCH] revert unrelated changes for the MR --- ...haracterize_Darks_NewDAQ_FastCCD_NBC.ipynb | 581 ++++++++ ...orrectionNotebook_NewDAQ_FastCCD_NBC.ipynb | 1279 +++++++++++++++++ .../FastCCD/PlotFromCalDB_FastCCD_NBC.ipynb | 505 +++++++ 3 files changed, 2365 insertions(+) create mode 100644 notebooks/FastCCD/Characterize_Darks_NewDAQ_FastCCD_NBC.ipynb create mode 100644 notebooks/FastCCD/CorrectionNotebook_NewDAQ_FastCCD_NBC.ipynb create mode 100644 notebooks/FastCCD/PlotFromCalDB_FastCCD_NBC.ipynb diff --git a/notebooks/FastCCD/Characterize_Darks_NewDAQ_FastCCD_NBC.ipynb b/notebooks/FastCCD/Characterize_Darks_NewDAQ_FastCCD_NBC.ipynb new file mode 100644 index 000000000..2f89983fe --- /dev/null +++ b/notebooks/FastCCD/Characterize_Darks_NewDAQ_FastCCD_NBC.ipynb @@ -0,0 +1,581 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# FastCCD Dark Characterization\n", + "\n", + "Author: I. KlaÄková, S. Hauf, Version 1.0\n", + "\n", + "The following notebook provides dark image analysis of the FastCCD detector.\n", + "\n", + "Dark characterization evaluates offset and noise of the detector and gives information about bad pixels. Resulting maps are saved as .h5 files for a latter use." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T10:54:38.999974Z", + "start_time": "2018-12-06T10:54:38.983406Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "in_folder = \"/gpfs/exfel/exp/SCS/201930/p900074/raw/\" # input folder, required\n", + "out_folder = 'gpfs/exfel/data/scratch/haufs/test/' # output folder, required\n", + "path_template = 'RAW-R{:04d}-DA05-S{{:05d}}.h5' # the template to use to access data\n", + "run = 321 # which run to read data from, required\n", + "number_dark_frames = 0 # number of images to be used, if set to 0 all available images are used\n", + "cluster_profile = \"noDB\" # ipcluster profile to use\n", + "operation_mode = \"FF\" #o r \"FF\". FS stands for frame-store and FF for full-frame opeartion\n", + "sigma_noise = 10. # Pixel exceeding 'sigmaNoise' * noise value in that pixel will be masked\n", + "h5path = '/INSTRUMENT/SCS_CDIDET_FCCD2M/DAQ/FCCD:daqOutput/data/image/pixels' # path in the HDF5 file the data is at\n", + "h5path_t = '/CONTROL/SCS_CDIDET_FCCD2M/CTRL/LSLAN/inputA/crdg/value' # path to find temperature at\n", + "h5path_cntrl = '/RUN/SCS_CDIDET_FCCD2M/DET/FCCD' # path to control data\n", + "cal_db_interface = \"tcp://max-exfl016:8020\" # calibration DB interface to use\n", + "local_output = False # output also in as H5 files\n", + "temp_limits = 5 # limits within which temperature is considered the same\n", + "sequence = 0 # sequence file to use\n", + "multi_iteration = False # use multiple iterations\n", + "use_dir_creation_date = True # use dir creation date\n", + "bad_pixel_offset_sigma = 5. # offset standard deviations above which to consider pixel bad \n", + "bad_pixel_noise_sigma = 5. # noise standard deviations above which to consider pixel bad \n", + "fix_temperature = 0. # fix temperature to this value, set to 0 to use slow control value" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T10:54:39.190907Z", + "start_time": "2018-12-06T10:54:39.186154Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions\n", + "from iCalibrationDB.detectors import DetectorTypes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T10:54:39.467334Z", + "start_time": "2018-12-06T10:54:39.427784Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "import XFELDetAna.xfelprofiler as xprof\n", + "\n", + "profiler = xprof.Profiler()\n", + "profiler.disable()\n", + "from XFELDetAna.util import env\n", + "env.iprofile = cluster_profile\n", + "\n", + "import warnings\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "from XFELDetAna import xfelpycaltools as xcal\n", + "from XFELDetAna import xfelpyanatools as xana\n", + "from XFELDetAna.plotting.util import prettyPlotting\n", + "prettyPlotting=True\n", + "from XFELDetAna.xfelreaders import ChunkReader\n", + "from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5\n", + "from cal_tools.tools import get_dir_creation_date\n", + "\n", + "import numpy as np\n", + "import h5py\n", + "import matplotlib.pyplot as plt\n", + "from iminuit import Minuit\n", + "\n", + "import time\n", + "import copy\n", + "\n", + "from prettytable import PrettyTable\n", + "\n", + "%matplotlib inline\n", + "\n", + "def nImagesOrLimit(nImages, limit):\n", + " if limit == 0:\n", + " return nImages\n", + " else:\n", + " return min(nImages, limit)\n", + " \n", + "sigmaNoise = sigma_noise" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T10:54:40.058101Z", + "start_time": "2018-12-06T10:54:40.042615Z" + }, + "collapsed": false + }, + "outputs": [], + "source": [ + "if operation_mode == \"FS\":\n", + " x = 960 # rows of the FastCCD to analyze in FS mode \n", + " y = 960 # columns of the FastCCD to analyze in FS mode \n", + " print('\\nYou are analyzing data in FS mode.')\n", + "else:\n", + " x = 1934 # rows of the FastCCD to analyze in FF mode \n", + " y = 960 # columns of the FastCCD to analyze in FF mode\n", + " print('\\nYou are analyzing data in FF mode.\\n')\n", + " \n", + "ped_dir = \"{}/r{:04d}\".format(in_folder, run)\n", + "fp_name = path_template.format(run)\n", + "\n", + "import datetime\n", + "creation_time = None\n", + "if use_dir_creation_date:\n", + " creation_time = get_dir_creation_date(in_folder, run)\n", + "\n", + "fp_path = '{}/{}'.format(ped_dir, fp_name)\n", + "\n", + "print(\"Reading data from: {}\\n\".format(fp_path))\n", + "print(\"Run is: {}\".format(run))\n", + "print(\"HDF5 path: {}\".format(h5path))\n", + "if creation_time:\n", + " print(\"Using {} as creation time\".format(creation_time.isoformat()))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T10:54:40.555804Z", + "start_time": "2018-12-06T10:54:40.452978Z" + }, + "collapsed": false + }, + "outputs": [], + "source": [ + "filename = fp_path.format(sequence)\n", + "sensorSize = [x, y]\n", + "chunkSize = 100 #Number of images to read per chunk\n", + "#Sensor area will be analysed according to blocksize\n", + "blockSize = [sensorSize[0]//2, sensorSize[1]//4] \n", + "xcal.defaultBlockSize = blockSize\n", + "cpuCores = 8 #Specifies the number of running cpu cores\n", + "memoryCells = 1 #FastCCD has 1 memory cell\n", + "#Specifies total number of images to proceed\n", + "nImages = fastccdreaderh5.getDataSize(filename, h5path)[0] \n", + "nImages = nImagesOrLimit(nImages, number_dark_frames)\n", + "print(\"\\nNumber of dark images to analyze: \",nImages)\n", + "commonModeBlockSize = blockSize\n", + "commonModeAxisR = 'row'#Axis along which common mode will be calculated\n", + "run_parallel = True\n", + "profile = False\n", + "\n", + "with h5py.File(filename, 'r') as f:\n", + " bias_voltage = int(f['{}/biasclock/bias/value'.format(h5path_cntrl)][0])\n", + " det_gain = int(f['{}/exposure/gain/value'.format(h5path_cntrl)][0])\n", + " integration_time = int(f['{}/acquisitionTime/value'.format(h5path_cntrl)][0])\n", + " temperature = np.mean(f[h5path_t])\n", + " temperature_k = temperature + 273.15\n", + " \n", + " if fix_temperature != 0.:\n", + " temperature_k = fix_temperature\n", + " print(\"Using fixed temperature\")\n", + " print(\"Bias voltage is {} V\".format(bias_voltage))\n", + " print(\"Detector gain is set to x{}\".format(det_gain))\n", + " print(\"Detector integration time is set to {}\".format(integration_time))\n", + " print(\"Mean temperature was {:0.2f} °C / {:0.2f} K\".format(temperature, temperature_k))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T10:54:41.584031Z", + "start_time": "2018-12-06T10:54:41.578462Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "reader = ChunkReader(filename, fastccdreaderh5.readData, \n", + " nImages, chunkSize, \n", + " path = h5path, \n", + " pixels_x = sensorSize[0],\n", + " pixels_y = sensorSize[1],)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T10:54:41.899511Z", + "start_time": "2018-12-06T10:54:41.864816Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "noiseCal = xcal.NoiseCalculator(sensorSize, memoryCells, \n", + " cores=cpuCores, blockSize=blockSize,\n", + " runParallel=run_parallel)\n", + "histCalRaw = xcal.HistogramCalculator(sensorSize, bins=1000, \n", + " range=[0, 10000], parallel=False, \n", + " memoryCells=memoryCells, \n", + " cores=cpuCores, blockSize=blockSize)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### First Iteration" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Characterization of dark images with purpose to create dark maps (offset, noise and bad pixel maps) is an iterative process. Firstly, initial offset and noise maps are produced from raw dark data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T10:55:21.238009Z", + "start_time": "2018-12-06T10:54:54.586435Z" + }, + "collapsed": false + }, + "outputs": [], + "source": [ + "for data in reader.readChunks():\n", + " data = np.bitwise_and(data.astype(np.uint16), 0b0011111111111111).astype(np.float32)\n", + " dx = np.count_nonzero(data, axis=(0, 1))\n", + " data = data[:,:,dx != 0]\n", + " histCalRaw.fill(data)\n", + " #Filling calculators with data\n", + " noiseCal.fill(data)\n", + " \n", + "offsetMap = noiseCal.getOffset() #Produce offset map\n", + "noiseMap = noiseCal.get() #Produce noise map\n", + "noiseCal.reset() #Reset noise calculator\n", + "print(\"Initial maps were created\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T10:56:20.686534Z", + "start_time": "2018-12-06T10:56:11.721829Z" + }, + "collapsed": false + }, + "outputs": [], + "source": [ + "#**************OFFSET MAP HISTOGRAM***********#\n", + "ho,co = np.histogram(offsetMap.flatten(), bins=700)\n", + "\n", + "do = {'x': co[:-1],\n", + " 'y': ho,\n", + " 'y_err': np.sqrt(ho[:]),\n", + " 'drawstyle': 'bars',\n", + " 'color': 'cornflowerblue',\n", + " }\n", + "\n", + "fig = xana.simplePlot(do, figsize='1col', aspect=2, \n", + " x_label = 'Offset (ADU)', \n", + " y_label=\"Counts\", y_log=True,\n", + " )\n", + " \n", + "\n", + "#*****NOISE MAP HISTOGRAM FROM THE OFFSET CORRECTED DATA*******#\n", + "hn,cn = np.histogram(noiseMap.flatten(), bins=200)\n", + "\n", + "dn = {'x': cn[:-1],\n", + " 'y': hn,\n", + " 'y_err': np.sqrt(hn[:]),\n", + " 'drawstyle': 'bars',\n", + " 'color': 'cornflowerblue',\n", + " }\n", + "\n", + "fig = xana.simplePlot(dn, figsize='1col', aspect=2, \n", + " x_label = 'Noise (ADU)', \n", + " y_label=\"Counts\", \n", + " y_log=True)\n", + "\n", + "\n", + "#**************HEAT MAPS*******************#\n", + "fig = xana.heatmapPlot(offsetMap[:,:,0],\n", + " x_label='Columns', y_label='Rows',\n", + " lut_label='Offset (ADU)',\n", + " x_range=(0,y),\n", + " y_range=(0,x), vmin=3000, vmax=4500)\n", + "\n", + "fig = xana.heatmapPlot(noiseMap[:,:,0],\n", + " x_label='Columns', y_label='Rows',\n", + " lut_label='Noise (ADU)',\n", + " x_range=(0,y),\n", + " y_range=(0,x), vmax=2*np.mean(noiseMap))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T10:56:22.741284Z", + "start_time": "2018-12-06T10:56:20.688393Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "\n", + "## offset\n", + "\n", + "metadata = ConstantMetaData()\n", + "offset = Constants.CCD(DetectorTypes.fastCCD).Offset()\n", + "offset.data = offsetMap.data\n", + "metadata.calibration_constant = offset\n", + "\n", + "# set the operating condition\n", + "condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,\n", + " integration_time=integration_time,\n", + " gain_setting=det_gain,\n", + " temperature=temperature_k,\n", + " pixels_x=1934,\n", + " pixels_y=960)\n", + "for parm in condition.parameters:\n", + " if parm.name == \"Sensor Temperature\":\n", + " parm.lower_deviation = temp_limits\n", + " parm.upper_deviation = temp_limits\n", + "\n", + "device = Detectors.fastCCD1\n", + "\n", + "\n", + "metadata.detector_condition = condition\n", + "\n", + "# specify the 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)\n", + "\n", + "## noise\n", + "\n", + "metadata = ConstantMetaData()\n", + "noise = Constants.CCD(DetectorTypes.fastCCD).Noise()\n", + "noise.data = noiseMap.data\n", + "metadata.calibration_constant = noise\n", + "\n", + "# set the operating condition\n", + "condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,\n", + " integration_time=integration_time,\n", + " gain_setting=det_gain,\n", + " temperature=temperature_k,\n", + " pixels_x=1934,\n", + " pixels_y=960)\n", + "\n", + "for parm in condition.parameters:\n", + " if parm.name == \"Sensor Temperature\":\n", + " parm.lower_deviation = temp_limits\n", + " parm.upper_deviation = temp_limits\n", + "\n", + "\n", + "device = Detectors.fastCCD1\n", + "\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)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "from cal_tools.enums import BadPixels\n", + "bad_pixels = np.zeros(offsetMap.shape, np.uint32)\n", + "mnoffset = np.nanmedian(offsetMap)\n", + "stdoffset = np.nanstd(offsetMap)\n", + "bad_pixels[(offsetMap < mnoffset-bad_pixel_offset_sigma*stdoffset) | \n", + " (offsetMap > mnoffset+bad_pixel_offset_sigma*stdoffset)] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value\n", + "\n", + "mnnoise = np.nanmedian(noiseMap)\n", + "stdnoise = np.nanstd(noiseMap)\n", + "bad_pixels[(noiseMap < mnnoise-bad_pixel_noise_sigma*stdnoise) | \n", + " (noiseMap > mnnoise+bad_pixel_noise_sigma*stdnoise)] = BadPixels.NOISE_OUT_OF_THRESHOLD.value\n", + "\n", + "fig = xana.heatmapPlot(np.log2(bad_pixels[:,:,0]),\n", + " x_label='Columns', y_label='Rows',\n", + " lut_label='Bad Pixel Value (ADU)',\n", + " x_range=(0,y),\n", + " y_range=(0,x), vmin=0, vmax=32)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "metadata = ConstantMetaData()\n", + "badpix = Constants.CCD(DetectorTypes.fastCCD).BadPixelsDark()\n", + "badpix.data = bad_pixels.data\n", + "metadata.calibration_constant = badpix\n", + "\n", + "# set the operating condition\n", + "condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,\n", + " integration_time=integration_time,\n", + " gain_setting=det_gain,\n", + " temperature=temperature_k,\n", + " pixels_x=1934,\n", + " pixels_y=960)\n", + "\n", + "for parm in condition.parameters:\n", + " if parm.name == \"Sensor Temperature\":\n", + " parm.lower_deviation = temp_limits\n", + " parm.upper_deviation = temp_limits\n", + "\n", + "\n", + "device = Detectors.fastCCD1\n", + "\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)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "histCalCorr = xcal.HistogramCalculator(sensorSize, bins=200, \n", + " range=[-200, 200], parallel=False, \n", + " memoryCells=memoryCells, \n", + " cores=cpuCores, blockSize=blockSize)\n", + "\n", + "\n", + "for data in reader.readChunks():\n", + " data = np.bitwise_and(data.astype(np.uint16), 0b0011111111111111).astype(np.float32)\n", + " data -= offsetMap.data\n", + " histCalCorr.fill(data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ho,eo,co,so = histCalCorr.get()\n", + "\n", + "\n", + "d = [{'x': co,\n", + " 'y': ho,\n", + " 'y_err': np.sqrt(ho[:]),\n", + " 'drawstyle': 'steps-mid',\n", + " 'errorstyle': 'bars',\n", + " 'errorcoarsing': 2,\n", + " 'label': 'Offset corr.'\n", + " },\n", + " \n", + " ]\n", + " \n", + "\n", + "fig = xana.simplePlot(d, aspect=1, x_label='Energy(ADU)', \n", + " y_label='Number of occurrences', figsize='2col',\n", + " y_log=True, x_range=(-50,500),\n", + " legend='top-center-frame-2col')" + ] + }, + { + "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" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/notebooks/FastCCD/CorrectionNotebook_NewDAQ_FastCCD_NBC.ipynb b/notebooks/FastCCD/CorrectionNotebook_NewDAQ_FastCCD_NBC.ipynb new file mode 100644 index 000000000..c80d17a47 --- /dev/null +++ b/notebooks/FastCCD/CorrectionNotebook_NewDAQ_FastCCD_NBC.ipynb @@ -0,0 +1,1279 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# FastCCD Data Correction ##\n", + "\n", + "Authors: I. KlaÄková, S. Hauf, Version 1.0\n", + "\n", + "The following notebook provides correction of images acquired with the FastCCD." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T15:54:23.218849Z", + "start_time": "2018-12-06T15:54:23.166497Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "in_folder = \"/gpfs/exfel/exp/SCS/201802/p002170/raw/\" # input folder, required\n", + "out_folder = '/gpfs/exfel/data/scratch/xcal/test/' # output folder, required\n", + "path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # path template in hdf5 file\n", + "path_inset = 'DA05'\n", + "run = 277 # run number\n", + "h5path = '/INSTRUMENT/SCS_CDIDET_FCCD2M/DAQ/FCCD:daqOutput/data/image' # path in HDF5 file\n", + "h5path_t = '/CONTROL/SCS_CDIDET_FCCD2M/CTRL/LSLAN/inputA/crdg/value' # temperature path in HDF5 file\n", + "h5path_cntrl = '/RUN/SCS_CDIDET_FCCD2M/DET/FCCD' # path to control data\n", + "cluster_profile = \"noDB\" #ipcluster profile to use\n", + "cpuCores = 16 #Specifies the number of running cpu cores\n", + "operation_mode = \"FF\" # FS stands for frame-store and FF for full-frame opeartion\n", + "split_evt_primary_threshold = 7. # primary threshold for split event classification in terms of n sigma noise\n", + "split_evt_secondary_threshold = 4. # secondary threshold for split event classification in terms of n sigma noise\n", + "split_evt_mip_threshold = 1000. # MIP threshold for event classification\n", + "cal_db_interface = \"tcp://max-exfl016:8015#8025\" # calibration DB interface to use\n", + "cal_db_timeout = 300000000 # timeout on caldb requests\n", + "sequences = [-1] # sequences to correct, set to -1 for all, range allowed\n", + "chunk_size_idim = 1 # H5 chunking size of output data\n", + "overwrite = True # overwrite existing files\n", + "do_pattern_classification = True # classify split events\n", + "sequences_per_node = 1 # sequences to correct per node\n", + "limit_images = 0 # limit images per file \n", + "correct_offset_drift = False # correct for offset drifts\n", + "use_dir_creation_date = True # use dir creation data for calDB queries\n", + "time_offset_days = 0 # offset in days for calibration parameters\n", + "photon_energy_gain_map = 2. # energy in keV\n", + "fix_temperature = 0. # fix temperature to this value, set to 0 to use slow control value\n", + "flipped_between = [\"2019-02-01\", \"2019-04-02\"] # detector was flipped during this timespan\n", + "temp_limits = 5 # limits within which temperature is considered the same\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, path_inset))\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", + " nsplits = len(seq_nums)//sequences_per_node+1\n", + " while nsplits > 8:\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)]\n", + " else:\n", + " return sequences" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T15:54:23.455376Z", + "start_time": "2018-12-06T15:54:23.413579Z" + } + }, + "outputs": [], + "source": [ + "import XFELDetAna.xfelprofiler as xprof\n", + "\n", + "profiler = xprof.Profiler()\n", + "profiler.disable()\n", + "from XFELDetAna.util import env\n", + "env.iprofile = cluster_profile\n", + "\n", + "import warnings\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "from XFELDetAna import xfelpycaltools as xcal\n", + "from XFELDetAna import xfelpyanatools as xana\n", + "from XFELDetAna.plotting.util import prettyPlotting\n", + "prettyPlotting=True\n", + "from XFELDetAna.xfelreaders import ChunkReader\n", + "from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5\n", + "\n", + "import numpy as np\n", + "import h5py\n", + "import matplotlib.pyplot as plt\n", + "from iminuit import Minuit\n", + "\n", + "import time\n", + "import copy\n", + "import os\n", + "\n", + "from prettytable import PrettyTable\n", + "\n", + "from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions\n", + "from iCalibrationDB.detectors import DetectorTypes\n", + "from cal_tools.tools import get_dir_creation_date\n", + "\n", + "from datetime import timedelta\n", + "\n", + "%matplotlib inline\n", + "\n", + "if sequences[0] == -1:\n", + " sequences = None\n", + " \n", + "offset_correction_args = (0.2459991787617141, 243.21639920846485)\n", + "t_base = 247.82\n", + "\n", + "if \"#\" in cal_db_interface:\n", + " prot, serv, ran = cal_db_interface.split(\":\")\n", + " r1, r2 = ran.split(\"#\")\n", + " cal_db_interface = \":\".join(\n", + " [prot, serv, str(np.random.randint(int(r1), int(r2)))])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T15:54:23.679069Z", + "start_time": "2018-12-06T15:54:23.662821Z" + } + }, + "outputs": [], + "source": [ + "if operation_mode == \"FS\":\n", + " x = 960 # rows of the FastCCD to analyze in FS mode \n", + " y = 960 # columns of the FastCCD to analyze in FS mode \n", + " print('\\nYou are analyzing data in FS mode.')\n", + "else:\n", + " x = 1934 # rows of the FastCCD to analyze in FF mode \n", + " y = 960 # columns of the FastCCD to analyze in FF mode\n", + " print('\\nYou are analyzing data in FF mode.')\n", + " \n", + "ped_dir = \"{}/r{:04d}\".format(in_folder, run)\n", + "out_folder = \"{}/r{:04d}\".format(out_folder, run)\n", + "fp_name = path_template.format(run, path_inset)\n", + "fp_path = '{}/{}'.format(ped_dir, fp_name)\n", + "\n", + "print(\"Reading data from: {}\\n\".format(fp_path))\n", + "print(\"Run is: {}\".format(run))\n", + "print(\"HDF5 path: {}\".format(h5path))\n", + "print(\"Data is output to: {}\".format(out_folder))\n", + "\n", + "import datetime\n", + "creation_time = None\n", + "if use_dir_creation_date:\n", + " creation_time = get_dir_creation_date(in_folder, run) + timedelta(days=time_offset_days)\n", + "if creation_time:\n", + " print(\"Using {} as creation time\".format(creation_time.isoformat()))\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T15:54:23.913269Z", + "start_time": "2018-12-06T15:54:23.868910Z" + } + }, + "outputs": [], + "source": [ + "\n", + "sensorSize = [x, y]\n", + "chunkSize = 100 #Number of images to read per chunk\n", + "blockSize = [sensorSize[0]//2, sensorSize[1]//4] #Sensor area will be analysed according to blocksize\n", + "xcal.defaultBlockSize = blockSize\n", + "memoryCells = 1 #FastCCD has 1 memory cell\n", + "#Specifies total number of images to proceed\n", + "\n", + "commonModeBlockSize = blockSize\n", + "commonModeAxisR = 'row'#Axis along which common mode will be calculated\n", + "run_parallel = True\n", + "profile = False\n", + "\n", + "temperature_k = 291\n", + "filename = fp_path.format(sequences[0] if sequences else 0)\n", + "with h5py.File(filename, 'r') as f:\n", + " bias_voltage = int(f['{}/biasclock/bias/value'.format(h5path_cntrl)][0])\n", + " det_gain = int(f['{}/exposure/gain/value'.format(h5path_cntrl)][0])\n", + " integration_time = int(f['{}/acquisitionTime/value'.format(h5path_cntrl)][0])\n", + " print(\"Bias voltage is {} V\".format(bias_voltage))\n", + " print(\"Detector gain is set to x{}\".format(det_gain))\n", + " print(\"Detector integration time is set to {}\".format(integration_time))\n", + " temperature = np.mean(f[h5path_t])\n", + " temperature_k = temperature + 273.15\n", + " if fix_temperature != 0.:\n", + " temperature_k = fix_temperature\n", + " print(\"Using fixed temperature\")\n", + " print(\"Mean temperature was {:0.2f} °C / {:0.2f} K at beginning of run\".format(temperature, temperature_k))\n", + " \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" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T15:54:24.088948Z", + "start_time": "2018-12-06T15:54:24.059925Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "dirlist = sorted(os.listdir(ped_dir))\n", + "file_list = []\n", + "total_sequences = 0\n", + "fsequences = []\n", + "for entry in dirlist:\n", + "\n", + " #only h5 file\n", + " abs_entry = \"{}/{}\".format(ped_dir, entry)\n", + " if os.path.isfile(abs_entry) and os.path.splitext(abs_entry)[1] == \".h5\":\n", + " \n", + " if sequences is None:\n", + " for seq in range(len(dirlist)):\n", + " \n", + " if path_template.format(run, path_inset).format(seq) in abs_entry:\n", + " file_list.append(abs_entry)\n", + " total_sequences += 1\n", + " fsequences.append(seq)\n", + " else:\n", + " for seq in sequences:\n", + " \n", + " if path_template.format(run, path_inset).format(seq) in abs_entry:\n", + " file_list.append(os.path.abspath(abs_entry))\n", + " total_sequences += 1\n", + " fsequences.append(seq)\n", + "sequences = fsequences" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T18:43:39.776018Z", + "start_time": "2018-12-06T18:43:39.759185Z" + } + }, + "outputs": [], + "source": [ + "import copy\n", + "from IPython.display import HTML, display, Markdown, Latex\n", + "import tabulate\n", + "print(\"Processing a total of {} sequence files\".format(total_sequences))\n", + "table = []\n", + "\n", + "\n", + "for k, f in enumerate(file_list):\n", + " table.append((k, f))\n", + "if len(table): \n", + " md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=[\"#\", \"file\"]))) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As a first step, dark maps have to be loaded." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T15:54:28.254544Z", + "start_time": "2018-12-06T15:54:24.709521Z" + } + }, + "outputs": [], + "source": [ + "offsetMap = None\n", + "badPixelMap = None\n", + "noiseMap = None\n", + "for i, g in enumerate([8, 2, 1]):\n", + " ## offset\n", + " metadata = ConstantMetaData()\n", + " offset = Constants.CCD(DetectorTypes.fastCCD).Offset()\n", + " metadata.calibration_constant = offset\n", + "\n", + " # set the operating condition\n", + " condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,\n", + " integration_time=integration_time,\n", + " gain_setting=g,\n", + " temperature=temperature_k,\n", + " pixels_x=1934,\n", + " pixels_y=960)\n", + "\n", + " for parm in condition.parameters:\n", + " if parm.name == \"Sensor Temperature\":\n", + " parm.lower_deviation = temp_limits\n", + " parm.upper_deviation = temp_limits\n", + "\n", + "\n", + " device = Detectors.fastCCD1\n", + "\n", + "\n", + " metadata.detector_condition = condition\n", + "\n", + "\n", + "\n", + " # specify the version for this constant\n", + " if creation_time is None:\n", + " metadata.calibration_constant_version = Versions.Now(device=device)\n", + " metadata.retrieve(cal_db_interface)\n", + " else:\n", + " metadata.calibration_constant_version = Versions.Timespan(device=device,\n", + " start=creation_time)\n", + " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)\n", + "\n", + "\n", + " if offsetMap is None:\n", + " offsetMap = np.zeros(list(offset.data.shape)+[3], np.float32)\n", + " offsetMap[...,i] = offset.data\n", + "\n", + " offset_temperature = None\n", + " for parm in condition.parameters:\n", + "\n", + " if parm.name == \"Sensor Temperature\":\n", + " offset_temperature = parm.value\n", + "\n", + " print(\"Temperature of detector when dark images (gain {}) for offset calculation \".format(g) +\n", + " \"were taken at: {:0.2f} K @ {}\".format(offset_temperature,\n", + " metadata.calibration_constant_version.begin_at))\n", + "\n", + " ## noise\n", + " metadata = ConstantMetaData()\n", + " noise = Constants.CCD(DetectorTypes.fastCCD).Noise()\n", + " metadata.calibration_constant = noise\n", + "\n", + " # set the operating condition\n", + " condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,\n", + " integration_time=integration_time,\n", + " gain_setting=g,\n", + " temperature=temperature_k,\n", + " pixels_x=1934,\n", + " pixels_y=960)\n", + "\n", + "\n", + " for parm in condition.parameters:\n", + " if parm.name == \"Sensor Temperature\":\n", + " parm.lower_deviation = temp_limits\n", + " parm.upper_deviation = temp_limits\n", + "\n", + "\n", + " device = Detectors.fastCCD1\n", + "\n", + "\n", + " metadata.detector_condition = condition\n", + "\n", + " # specify the version for this constant\n", + " if creation_time is None:\n", + " metadata.calibration_constant_version = Versions.Now(device=device)\n", + " metadata.retrieve(cal_db_interface)\n", + " else:\n", + " metadata.calibration_constant_version = Versions.Timespan(device=device,\n", + " start=creation_time)\n", + " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)\n", + "\n", + " if noiseMap is None:\n", + " noiseMap = np.zeros(list(noise.data.shape)+[3], np.float32)\n", + " noiseMap[...,i] = noise.data\n", + "\n", + "\n", + " ## bad pixels \n", + "\n", + " metadata = ConstantMetaData()\n", + " bpix = Constants.CCD(DetectorTypes.fastCCD).BadPixelsDark()\n", + " metadata.calibration_constant = bpix\n", + "\n", + " # set the operating condition\n", + " condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,\n", + " integration_time=integration_time,\n", + " gain_setting=g,\n", + " temperature=temperature_k,\n", + " pixels_x=1934,\n", + " pixels_y=960)\n", + "\n", + " for parm in condition.parameters:\n", + " if parm.name == \"Sensor Temperature\":\n", + " parm.lower_deviation = temp_limits\n", + " parm.upper_deviation = temp_limits\n", + "\n", + "\n", + " device = Detectors.fastCCD1\n", + "\n", + "\n", + " metadata.detector_condition = condition\n", + "\n", + " # specify the version for this constant\n", + " if creation_time is None:\n", + " metadata.calibration_constant_version = Versions.Now(device=device)\n", + " metadata.retrieve(cal_db_interface)\n", + " else:\n", + " metadata.calibration_constant_version = Versions.Timespan(device=device,\n", + " start=creation_time)\n", + " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)\n", + "\n", + " if badPixelMap is None:\n", + " badPixelMap = np.zeros(list(bpix.data.shape)+[3], np.uint32)\n", + " badPixelMap[...,i] = bpix.data\n", + " \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Loading cti and relative gain values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T15:54:28.343869Z", + "start_time": "2018-12-06T15:54:28.271344Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "## relative gain\n", + "\n", + "metadata = ConstantMetaData()\n", + "relgain = Constants.CCD(DetectorTypes.fastCCD).RelativeGain()\n", + "metadata.calibration_constant = relgain\n", + "\n", + "# set the operating condition\n", + "condition = Conditions.Illuminated.CCD(bias_voltage=bias_voltage,\n", + " integration_time=integration_time,\n", + " gain_setting=0,\n", + " temperature=temperature_k,\n", + " pixels_x=1934,\n", + " pixels_y=960, photon_energy=photon_energy_gain_map)\n", + "device = Detectors.fastCCD1\n", + "\n", + "\n", + "metadata.detector_condition = condition\n", + "\n", + "# specify the a version for this constant\n", + "metadata.calibration_constant_version = Versions.Now(device=device)\n", + "metadata.retrieve(cal_db_interface)\n", + "\n", + "relGain = relgain.data[::-1,...]\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "relGainCA = copy.copy(relGain)\n", + "relGainC = relGainCA[:relGainCA.shape[0]//2,...]\n", + "ctiA = np.ones(relGainCA.shape[:2])\n", + "cti = np.ones(relGainC.shape[:2])\n", + "i = 0\n", + "idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)\n", + "mn1 = np.nanmean(relGainC[i, ~idx, 0])\n", + "\n", + "for i in range(1, relGainC.shape[0]):\n", + " idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)\n", + " mn2 = np.nanmean(relGainC[i, ~idx, 0])\n", + " cti[i,:] = mn2/mn1\n", + "ctiA[:relGainCA.shape[0]//2,...] = cti\n", + "\n", + "relGainC = relGainCA[relGainCA.shape[0]//2:,...]\n", + "\n", + "\n", + "cti = np.ones(relGainC.shape[:2])\n", + "i = -1\n", + "idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)\n", + "mn1 = np.nanmean(relGainC[i, ~idx, 0])\n", + "\n", + "for i in range(relGainC.shape[0]-1, 1, -1):\n", + " idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)\n", + " mn2 = np.nanmean(relGainC[i, ~idx, 0])\n", + " cti[i,:] = mn2/mn1\n", + "\n", + "ctiA[relGainCA.shape[0]//2:,...] = cti\n", + "\n", + "relGainCA = copy.copy(relGain)\n", + "relGainC = relGainCA[:relGainCA.shape[0]//2,...]\n", + "for i in range(relGainC.shape[1]):\n", + " idx = (relGainC[:,i, 0] < 0.95) | (relGainC[:,i,0] > 1.05)\n", + " relGainC[idx,i,0] = np.nanmean(relGainC[~idx,i,0])\n", + " relGainC[idx,i,1] = np.nanmean(relGainC[~idx,i,1])\n", + " relGainC[idx,i,2] = np.nanmean(relGainC[~idx,i,2])\n", + "relGainCA[:relGainCA.shape[0]//2,...] = relGainC\n", + "relGainC = relGainCA[relGainCA.shape[0]//2:,...]\n", + "for i in range(relGainC.shape[1]):\n", + " idx = (relGainC[:,i, 0] < 0.95) | (relGainC[:,i,0] > 1.05)\n", + " relGainC[idx,i,0] = np.nanmean(relGainC[~idx,i,0])\n", + " relGainC[idx,i,1] = np.nanmean(relGainC[~idx,i,1])\n", + " relGainC[idx,i,2] = np.nanmean(relGainC[~idx,i,2])\n", + "relGainCA[relGainCA.shape[0]//2:,...] = relGainC\n", + "relGainC = relGainCA*ctiA[...,None]\n", + "\n", + "relGain = relGainC" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import dateutil.parser\n", + "flipped_between = [dateutil.parser.parse(d) for d in flipped_between]\n", + "flip_rgain = creation_time >= flipped_between[0] and creation_time <= flipped_between[1]\n", + "flip_rgain &= (metadata.calibration_constant_version.begin_at.replace(tzinfo=None) >= flipped_between[0] \n", + " and metadata.calibration_constant_version.begin_at.replace(tzinfo=None) <= flipped_between[1])\n", + "print(\"Accounting for flipped detector: {}\".format(flip_rgain))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T15:54:28.771629Z", + "start_time": "2018-12-06T15:54:28.346051Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "#************************Calculators************************#\n", + "\n", + "\n", + "cmCorrection = xcal.CommonModeCorrection([x, y], \n", + " commonModeBlockSize, \n", + " commonModeAxisR,\n", + " nCells = memoryCells, \n", + " noiseMap = noiseMap,\n", + " runParallel=True,\n", + " stats=True)\n", + "\n", + "patternClassifierLH = xcal.PatternClassifier([x//2, y], \n", + " noiseMap[:x//2, :], \n", + " split_evt_primary_threshold, \n", + " split_evt_secondary_threshold,\n", + " split_evt_mip_threshold,\n", + " tagFirstSingles = 0, \n", + " nCells=memoryCells, \n", + " cores=cpuCores, \n", + " allowElongated = False,\n", + " blockSize=[x//2, y],\n", + " runParallel=True)\n", + "\n", + "\n", + "\n", + "patternClassifierUH = xcal.PatternClassifier([x//2, y], \n", + " noiseMap[x//2:, :], \n", + " split_evt_primary_threshold, \n", + " split_evt_secondary_threshold,\n", + " split_evt_mip_threshold,\n", + " tagFirstSingles = 0, \n", + " nCells=memoryCells, \n", + " cores=cpuCores, \n", + " allowElongated = False,\n", + " blockSize=[x//2, y],\n", + " runParallel=True)\n", + "\n", + " \n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:08:51.886343Z", + "start_time": "2018-12-06T16:08:51.842837Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "#*****************Histogram Calculators******************#\n", + "\n", + "histCalOffsetCor = xcal.HistogramCalculator([x, y], \n", + " bins=500, \n", + " range=[-50, 1000],\n", + " nCells=memoryCells, \n", + " cores=cpuCores,\n", + " blockSize=blockSize)\n", + "\n", + "histCalPcorr = xcal.HistogramCalculator([x, y], \n", + " bins=500, \n", + " range=[-50, 1000],\n", + " nCells=memoryCells, \n", + " cores=cpuCores,\n", + " blockSize=blockSize)\n", + "\n", + "histCalPcorrS = xcal.HistogramCalculator([x, y], \n", + " bins=500, \n", + " range=[-50, 1000],\n", + " nCells=memoryCells, \n", + " cores=cpuCores,\n", + " blockSize=blockSize)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Applying corrections" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:08:52.441784Z", + "start_time": "2018-12-06T16:08:52.437284Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "patternClassifierLH._imagesPerChunk = 500\n", + "patternClassifierUH._imagesPerChunk = 500\n", + "patternClassifierLH.debug()\n", + "patternClassifierUH.debug()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:08:53.042555Z", + "start_time": "2018-12-06T16:08:53.034522Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "histCalOffsetCor.debug()\n", + "histCalPcorr.debug()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:08:53.551111Z", + "start_time": "2018-12-06T16:08:53.531064Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "def copy_and_sanitize_non_cal_data(infile, outfile, h5base):\n", + " \n", + " if h5base.startswith(\"/\"):\n", + " h5base = h5base[1:]\n", + " dont_copy = ['pixels']\n", + " dont_copy = [h5base+\"/{}\".format(do)\n", + " for do in dont_copy]\n", + "\n", + " def visitor(k, item):\n", + " if k not in dont_copy:\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)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:10:55.917179Z", + "start_time": "2018-12-06T16:09:01.603633Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "mean_im = None\n", + "single_im = None\n", + "mean_im_cc = None\n", + "single_im_cc = None\n", + "drift_lh = []\n", + "drift_uh = []\n", + "offsetMap = np.squeeze(offsetMap)\n", + "noiseMap = np.squeeze(noiseMap)\n", + "badPixelMap = np.squeeze(badPixelMap)\n", + "relGain = np.squeeze(relGain)\n", + "for k, f in enumerate(file_list):\n", + " with h5py.File(f, 'r', driver='core') as infile:\n", + " out_fileb = \"{}/{}\".format(out_folder, f.split(\"/\")[-1])\n", + " out_file = out_fileb.replace(\"RAW\", \"CORR\")\n", + " #out_filed = out_fileb.replace(\"RAW\", \"CORR-SC\")\n", + "\n", + " data = None\n", + " noise = None\n", + " try:\n", + " with h5py.File(out_file, \"w\") as ofile:\n", + " \n", + " copy_and_sanitize_non_cal_data(infile, ofile, h5path)\n", + " data = infile[h5path+\"/pixels\"][()]\n", + " nzidx = np.count_nonzero(data, axis=(1,2))\n", + " data = data[nzidx != 0, ...]\n", + " if limit_images > 0:\n", + " data = data[:limit_images,...]\n", + " oshape = data.shape\n", + " data = np.moveaxis(data, 0, 2)\n", + " ddset = ofile.create_dataset(h5path+\"/pixels\",\n", + " oshape,\n", + " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", + " dtype=np.float32)\n", + " \n", + " ddsetm = ofile.create_dataset(h5path+\"/mask\",\n", + " oshape,\n", + " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", + " dtype=np.uint32, compression=\"gzip\")\n", + " \n", + " ddsetg = ofile.create_dataset(h5path+\"/gain\",\n", + " oshape,\n", + " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", + " dtype=np.uint8, compression=\"gzip\")\n", + " \n", + " gain = np.right_shift(data, 14)\n", + " \n", + " gain[gain != 0] -= 1\n", + " \n", + " fstride = 1\n", + " if not flip_rgain: # rgain was taken during flipped orientation\n", + " fstride = -1\n", + " \n", + " data = np.bitwise_and(data, 0b0011111111111111).astype(np.float32) \n", + " omap = np.repeat(offsetMap[...,None,:], data.shape[2], axis=2)\n", + " rmap = np.repeat(relGain[:,::fstride,None,:], data.shape[2], axis=2)\n", + " nmap = np.repeat(noiseMap[...,None,:], data.shape[2], axis=2)\n", + " bmap = np.repeat(badPixelMap[...,None,:], data.shape[2], axis=2)\n", + " offset = np.choose(gain, (omap[...,0], omap[...,1], omap[...,2]))\n", + " rg = np.choose(gain, (rmap[...,0], rmap[...,1], rmap[...,2]))\n", + " noise = np.choose(gain, (nmap[...,0], nmap[...,1], nmap[...,2]))\n", + " bpix = np.choose(gain, (bmap[...,0], bmap[...,1], bmap[...,2]))\n", + " \n", + " data -= offset\n", + " data *= rg\n", + " \n", + " if correct_offset_drift:\n", + " lhd = np.mean(data[x//2-10:x//2,y//2-5:y//2+5,:], axis=(0,1))\n", + " data[:x//2, :, :] -= lhd\n", + " drift_lh.append(lhd)\n", + " \n", + " uhd = np.mean(data[x//2:x//2+10,y//2-5:y//2+5,:], axis=(0,1)) \n", + " data[x//2:, :, :] -= uhd\n", + " drift_uh.append(lhd)\n", + " \n", + " histCalOffsetCor.fill(data)\n", + "\n", + " \n", + " ddset[...] = np.moveaxis(data, 2, 0)\n", + " ddsetm[...] = np.moveaxis(bpix, 2, 0)\n", + " ddsetg[...] = np.moveaxis(gain, 2, 0).astype(np.uint8)\n", + " \n", + " if mean_im is None:\n", + " mean_im = np.nanmean(data, axis=2)\n", + " single_im = data[...,0]\n", + " \n", + " if do_pattern_classification:\n", + " \n", + " ddsetcm = ofile.create_dataset(h5path+\"/pixels_cm\",\n", + " oshape,\n", + " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", + " dtype=np.float32)\n", + "\n", + " ddsetc = ofile.create_dataset(h5path+\"/pixels_classified\",\n", + " oshape,\n", + " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", + " dtype=np.float32, compression=\"gzip\")\n", + "\n", + " ddsetp = ofile.create_dataset(h5path+\"/patterns\",\n", + " oshape,\n", + " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", + " dtype=np.int32, compression=\"gzip\")\n", + " \n", + "\n", + " patternClassifierLH._noisemap = noise[:x//2, :, :]\n", + " patternClassifierUH._noisemap = noise[x//2:, :, :]\n", + "\n", + " data = cmCorrection.correct(data) # correct for the row common mode\n", + " ddsetcm[...] = np.moveaxis(data, 2, 0)\n", + "\n", + " dataLH = data[:x//2, :, :]\n", + " dataUH = data[x//2:, :, :]\n", + "\n", + " dataLH, patternsLH = patternClassifierLH.classify(dataLH)\n", + " dataUH, patternsUH = patternClassifierUH.classify(dataUH)\n", + "\n", + " data[:x//2, :, :] = dataLH\n", + " data[x//2:, :, :] = dataUH\n", + "\n", + " patterns = np.zeros(data.shape, patternsLH.dtype)\n", + " patterns[:x//2, :, :] = patternsLH\n", + " patterns[x//2:, :, :] = patternsUH\n", + "\n", + " data[data < split_evt_primary_threshold*noise] = 0\n", + " ddsetc[...] = np.moveaxis(data, 2, 0)\n", + " ddsetp[...] = np.moveaxis(patterns, 2, 0)\n", + "\n", + " histCalPcorr.fill(data)\n", + " data[patterns != 100] = np.nan\n", + " histCalPcorrS.fill(data)\n", + "\n", + " if mean_im_cc is None:\n", + " mean_im_cc = np.nanmean(data, axis=2)\n", + " single_im_cc = data[...,0]\n", + " \n", + " except Exception as e:\n", + " print(\"Couldn't calibrate data in {}: {}\".format(f, e))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:10:56.094985Z", + "start_time": "2018-12-06T16:10:55.918900Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "if correct_offset_drift:\n", + " lhds = np.concatenate(drift_lh)\n", + " uhds = np.concatenate(drift_uh)\n", + " fig = plt.figure(figsize=(10,5))\n", + " ax = fig.add_subplot(111)\n", + " ax.plot(lhds, label=\"Lower hem.\")\n", + " ax.plot(uhds, label=\"Upper hem.\")\n", + " ax.set_xlabel(\"Frame #\")\n", + " ax.set_xlabel(\"Offset drift (ADU)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:10:56.126409Z", + "start_time": "2018-12-06T16:10:56.096242Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "if do_pattern_classification:\n", + " print(\"******************LOWER HEMISPHERE******************\\n\")\n", + "\n", + " patternStatsLH = patternClassifierLH.getPatternStats()\n", + " fig = plt.figure(figsize=(15,15))\n", + " ax = fig.add_subplot(4,4,1)\n", + " sfields = [\"singles\", \"first singles\", \"clusters\"]\n", + " mfields = [\"doubles\", \"triples\", \"quads\"]\n", + " relativeOccurances = []\n", + " labels = []\n", + " for i, f in enumerate(sfields):\n", + " relativeOccurances.append(patternStatsLH[f])\n", + " labels.append(f)\n", + " for i, f in enumerate(mfields):\n", + " for k in range(len(patternStatsLH[f])):\n", + " relativeOccurances.append(patternStatsLH[f][k])\n", + " labels.append(f+\"(\"+str(k)+\")\")\n", + " relativeOccurances = np.array(relativeOccurances, np.float)\n", + " relativeOccurances/=np.sum(relativeOccurances)\n", + " pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)\n", + " ax.set_title(\"Pattern occurrence\")\n", + " # Set aspect ratio to be equal so that pie is drawn as a circle.\n", + " a = ax.axis('equal')\n", + "\n", + " smaps = [\"singlemap\", \"firstsinglemap\", \"clustermap\"]\n", + " for i, m in enumerate(smaps):\n", + "\n", + " ax = fig.add_subplot(4,4,2+i)\n", + "\n", + " pmap = ax.imshow(patternStatsLH[m], interpolation=\"nearest\", vmax=2*np.nanmedian(patternStatsLH[m]))\n", + " ax.set_title(m)\n", + " cb = fig.colorbar(pmap)\n", + "\n", + " mmaps = [\"doublemap\", \"triplemap\", \"quadmap\"]\n", + " k = 0\n", + " for i, m in enumerate(mmaps):\n", + "\n", + " for j in range(4):\n", + " ax = fig.add_subplot(4,4,2+len(smaps)+k)\n", + " pmap = ax.imshow(patternStatsLH[m][j], interpolation=\"nearest\", vmax=2*np.median(patternStatsLH[m][j]))\n", + " ax.set_title(m+\"(\"+str(j)+\")\")\n", + " cb = fig.colorbar(pmap)\n", + " k+=1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:10:56.176160Z", + "start_time": "2018-12-06T16:10:56.127853Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "if do_pattern_classification:\n", + " patternStatsUH = patternClassifierUH.getPatternStats()\n", + " fig = plt.figure(figsize=(15,15))\n", + " ax = fig.add_subplot(4,4,1)\n", + " sfields = [\"singles\", \"first singles\", \"clusters\"]\n", + " mfields = [\"doubles\", \"triples\", \"quads\"]\n", + " relativeOccurances = []\n", + " labels = []\n", + " for i, f in enumerate(sfields):\n", + " relativeOccurances.append(patternStatsUH[f])\n", + " labels.append(f)\n", + " for i, f in enumerate(mfields):\n", + " for k in range(len(patternStatsUH[f])):\n", + " relativeOccurances.append(patternStatsUH[f][k])\n", + " labels.append(f+\"(\"+str(k)+\")\")\n", + " relativeOccurances = np.array(relativeOccurances, np.float)\n", + " relativeOccurances/=np.sum(relativeOccurances)\n", + " pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)\n", + " ax.set_title(\"Pattern occurrence\")\n", + " # Set aspect ratio to be equal so that pie is drawn as a circle.\n", + " a = ax.axis('equal')\n", + "\n", + " smaps = [\"singlemap\", \"firstsinglemap\", \"clustermap\"]\n", + " for i, m in enumerate(smaps):\n", + "\n", + " ax = fig.add_subplot(4,4,2+i)\n", + "\n", + " pmap = ax.imshow(patternStatsUH[m], interpolation=\"nearest\", vmax=2*np.nanmedian(patternStatsUH[m]))\n", + " ax.set_title(m)\n", + " cb = fig.colorbar(pmap)\n", + "\n", + " mmaps = [\"doublemap\", \"triplemap\", \"quadmap\"]\n", + " k = 0\n", + " for i, m in enumerate(mmaps):\n", + "\n", + " for j in range(4):\n", + " ax = fig.add_subplot(4,4,2+len(smaps)+k)\n", + " pmap = ax.imshow(patternStatsUH[m][j], interpolation=\"nearest\", vmax=np.median(patternStatsUH[m][j]))\n", + " ax.set_title(m+\"(\"+str(j)+\")\")\n", + " cb = fig.colorbar(pmap)\n", + " k+=1\n", + "\n", + " print(\"******************UPPER HEMISPHERE******************\\n\") " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:10:56.190150Z", + "start_time": "2018-12-06T16:10:56.177570Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "if do_pattern_classification:\n", + " t0 = PrettyTable()\n", + " t0.title = \"Total number of Counts after all corrections\"\n", + " t0.field_names = [\"Hemisphere\",\"Singles\", \"First-singles\", \"Clusters\"]\n", + " t0.add_row([\"LH\", patternStatsLH['singles'], patternStatsLH['first singles'], patternStatsLH['clusters']])\n", + " t0.add_row([\"UH\", patternStatsUH['singles'], patternStatsUH['first singles'], patternStatsUH['clusters']])\n", + "\n", + " print(t0)\n", + "\n", + " t1 = PrettyTable()\n", + "\n", + " t1.field_names = [\"Index\",\"D-LH\", \"D-UH\", \"T-LH\", \"T-UH\", \"Q-LH\", \"Q-UH\"]\n", + "\n", + " t1.add_row([0, patternStatsLH['doubles'][0], patternStatsUH['doubles'][0], patternStatsLH['triples'][0], patternStatsUH['triples'][0], patternStatsLH['quads'][0], patternStatsUH['quads'][0]])\n", + " t1.add_row([1, patternStatsLH['doubles'][1], patternStatsUH['doubles'][1], patternStatsLH['triples'][1], patternStatsUH['triples'][1], patternStatsLH['quads'][1], patternStatsUH['quads'][1]])\n", + " t1.add_row([2, patternStatsLH['doubles'][2], patternStatsUH['doubles'][2], patternStatsLH['triples'][2], patternStatsUH['triples'][2], patternStatsLH['quads'][2], patternStatsUH['quads'][2]])\n", + " t1.add_row([3, patternStatsLH['doubles'][3], patternStatsUH['doubles'][3], patternStatsLH['triples'][3], patternStatsUH['triples'][3], patternStatsLH['quads'][3], patternStatsUH['quads'][3]])\n", + "\n", + " print(t1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:10:56.203219Z", + "start_time": "2018-12-06T16:10:56.191509Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "if do_pattern_classification:\n", + " doublesLH = patternStatsLH['doubles'][0] + patternStatsLH['doubles'][1] + patternStatsLH['doubles'][2] + patternStatsLH['doubles'][3]\n", + " triplesLH = patternStatsLH['triples'][0] + patternStatsLH['triples'][1] + patternStatsLH['triples'][2] + patternStatsLH['triples'][3]\n", + " quadsLH = patternStatsLH['quads'][0] + patternStatsLH['quads'][1] + patternStatsLH['quads'][2] + patternStatsLH['quads'][3]\n", + " allsinglesLH = patternStatsLH['singles'] + patternStatsLH['first singles']\n", + " eventsLH = allsinglesLH + doublesLH + triplesLH + quadsLH\n", + "\n", + " doublesUH = patternStatsUH['doubles'][0] + patternStatsUH['doubles'][1] + patternStatsUH['doubles'][2] + patternStatsUH['doubles'][3]\n", + " triplesUH = patternStatsUH['triples'][0] + patternStatsUH['triples'][1] + patternStatsUH['triples'][2] + patternStatsUH['triples'][3]\n", + " quadsUH = patternStatsUH['quads'][0] + patternStatsUH['quads'][1] + patternStatsUH['quads'][2] + patternStatsUH['quads'][3]\n", + " allsinglesUH = patternStatsUH['singles'] + patternStatsUH['first singles']\n", + " eventsUH = allsinglesUH + doublesUH + triplesUH + quadsUH\n", + "\n", + " reloccurLH = np.array([allsinglesLH/eventsLH, doublesLH/eventsLH, triplesLH/eventsLH, quadsLH/eventsLH])\n", + " reloccurUH = np.array([allsinglesUH/eventsUH, doublesUH/eventsUH, triplesUH/eventsUH, quadsUH/eventsUH])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:10:56.212586Z", + "start_time": "2018-12-06T16:10:56.204731Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "if do_pattern_classification:\n", + " fig = plt.figure(figsize=(10,5))\n", + " ax = fig.add_subplot(1,2,1)\n", + " labels = ['singles', 'doubles', 'triples', 'quads']\n", + " pie = ax.pie(reloccurLH, labels=labels, autopct='%1.1f%%', shadow=True)\n", + " ax.set_title(\"Pattern occurrence LH\")\n", + " # Set aspect ratio to be equal so that pie is drawn as a circle.\n", + " a = ax.axis('equal')\n", + " ax = fig.add_subplot(1,2,2)\n", + " pie = ax.pie(reloccurUH, labels=labels, autopct='%1.1f%%', shadow=True)\n", + " ax.set_title(\"Pattern occurrence UH\")\n", + " # Set aspect ratio to be equal so that pie is drawn as a circle.\n", + " a = ax.axis('equal')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:13:12.889583Z", + "start_time": "2018-12-06T16:13:11.122653Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "ho,eo,co,so = histCalOffsetCor.get()\n", + "\n", + "\n", + "d = [{'x': co,\n", + " 'y': ho,\n", + " 'y_err': np.sqrt(ho[:]),\n", + " 'drawstyle': 'steps-mid',\n", + " 'errorstyle': 'bars',\n", + " 'errorcoarsing': 2,\n", + " 'label': 'Offset corr.'\n", + " },\n", + " \n", + " ]\n", + " \n", + "\n", + "fig = xana.simplePlot(d, aspect=1, x_label='Energy(ADU)', \n", + " y_label='Number of occurrences', figsize='2col',\n", + " y_log=True, x_range=(-50,500),\n", + " legend='top-center-frame-2col')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:12:57.289742Z", + "start_time": "2018-12-06T16:12:45.529734Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "if do_pattern_classification:\n", + " h1,e1L,c1L,s1L = histCalPcorr.get()\n", + " h1s,e1Ls,c1Ls,s1Ls = histCalPcorrS.get()\n", + "\n", + "\n", + " d = [\n", + " {'x': c1L,\n", + " 'y': h1,\n", + " 'y_err': np.sqrt(h1[:]),\n", + " 'drawstyle': 'steps-mid',\n", + " 'label': 'Split event corrected'},\n", + " {'x': c1Ls,\n", + " 'y': h1s,\n", + " 'y_err': np.sqrt(h1s[:]),\n", + " 'drawstyle': 'steps-mid',\n", + " 'label': 'Single pixel hits'}\n", + " ]\n", + "\n", + "\n", + " fig = xana.simplePlot(d, aspect=1, x_label='Energy(ADU)', \n", + " y_label='Number of occurrences', figsize='2col',\n", + " y_log=True, x_range=(0,200),x_log=False,\n", + " legend='top-center-frame-2col')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Mean Image of first Sequence ##" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:11:08.317130Z", + "start_time": "2018-12-06T16:11:05.788655Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "fig = xana.heatmapPlot(mean_im,\n", + " x_label='Columns', y_label='Rows',\n", + " lut_label='Signal (ADU)',\n", + " x_range=(0,y),\n", + " y_range=(0,x), vmin=-50, vmax=500)\n", + "\n", + "if do_pattern_classification:\n", + " fig = xana.heatmapPlot(mean_im_cc,\n", + " x_label='Columns', y_label='Rows',\n", + " lut_label='Signal (ADU)',\n", + " x_range=(0,y),\n", + " y_range=(0,x), vmin=-50, vmax=500)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "## Single Shot of first Sequnce ##" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "ExecuteTime": { + "end_time": "2018-12-06T16:11:10.908912Z", + "start_time": "2018-12-06T16:11:08.318486Z" + }, + "collapsed": true + }, + "outputs": [], + "source": [ + "fig = xana.heatmapPlot(single_im,\n", + " x_label='Columns', y_label='Rows',\n", + " lut_label='Signal (ADU)',\n", + " x_range=(0,y),\n", + " y_range=(0,x), vmin=-50, vmax=500)\n", + "\n", + "if do_pattern_classification:\n", + " fig = xana.heatmapPlot(single_im_cc,\n", + " x_label='Columns', y_label='Rows',\n", + " lut_label='Signal (ADU)',\n", + " x_range=(0,y),\n", + " y_range=(0,x), vmin=-50, vmax=500)" + ] + }, + { + "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.7" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/notebooks/FastCCD/PlotFromCalDB_FastCCD_NBC.ipynb b/notebooks/FastCCD/PlotFromCalDB_FastCCD_NBC.ipynb new file mode 100644 index 000000000..908d75fea --- /dev/null +++ b/notebooks/FastCCD/PlotFromCalDB_FastCCD_NBC.ipynb @@ -0,0 +1,505 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Statistical analysis of calibration factors#\n", + "\n", + "Author: Mikhail Karnevskiy, Steffen Hauf, Version 0.1\n", + "\n", + "A description of the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cluster_profile = \"noDB\" # The ipcluster profile to use\n", + "start_date = \"2019-01-30\" # date to start investigation interval from\n", + "end_date = \"2019-08-30\" # date to end investigation interval at, can be \"now\"\n", + "nconstants = 10 # Number of time stamps to plot. If not 0, overcome start_date.\n", + "dclass=\"CCD\" # Detector class\n", + "db_module = \"fastCCD1\" # detector entry in the DB to investigate\n", + "constants = [\"Noise\", \"Offset\"] # constants to plot\n", + "\n", + "gain_setting = [0,1,2,8] # gain stages\n", + "bias_voltage = [79] # Bias voltage\n", + "temperature = [235, 216, 245] # Operation temperature\n", + "integration_time = [1, 50] # Integration time\n", + "pixels_x=[1934] # number of pixels along X axis\n", + "pixels_y=[960] # number of pixels along Y axis\n", + "max_time = 15 # max time margin in minutes to match bad pixels\n", + "parameter_names = ['bias_voltage', 'integration_time', 'temperature', \n", + " 'gain_setting', 'pixels_x', 'pixels_y'] # names of parameters\n", + "\n", + "separate_plot = ['integration_time', 'gain_setting', 'temperature'] # Plot on separate plots\n", + "photon_energy = 9.2 # Photon energy of the beam\n", + "out_folder = \"/gpfs/exfel/data/scratch/karnem/test_FCCD/\" # output folder\n", + "use_existing = \"\" # If not empty, constants stored in given folder will be used\n", + "cal_db_interface = \"tcp://max-exfl016:8015#8025\" # the database interface to use\n", + "cal_db_timeout = 180000 # timeout on caldb requests\",\n", + "plot_range = 3 # range for plotting in units of median absolute deviations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "import copy\n", + "import datetime\n", + "import dateutil.parser\n", + "import numpy as np\n", + "from operator import itemgetter\n", + "import os\n", + "import sys\n", + "import warnings\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "from iCalibrationDB import Constants, Conditions, Detectors, ConstantMetaData\n", + "from cal_tools.tools import get_from_db\n", + "from cal_tools.ana_tools import (save_dict_to_hdf5, load_data_from_hdf5, \n", + " HMType, hm_combine,\n", + " combine_lists, get_range)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Prepare variables\n", + "spShape = (967, 10) # Shape of superpixel\n", + "\n", + "parameters = [globals()[x] for x in parameter_names]\n", + "\n", + "constantsDark = {'Noise': 'BadPixelsDark',\n", + " 'Offset': 'BadPixelsDark'}\n", + "print('Bad pixels data: ', constantsDark)\n", + "\n", + "# Define parameters in order to perform loop over time stamps\n", + "start = datetime.datetime.now() if start_date.upper() == \"NOW\" else dateutil.parser.parse(\n", + " start_date)\n", + "end = datetime.datetime.now() if end_date.upper() == \"NOW\" else dateutil.parser.parse(\n", + " end_date)\n", + "\n", + "# Create output folder\n", + "os.makedirs(out_folder, exist_ok=True)\n", + "\n", + "# Get getector conditions\n", + "det = getattr(Detectors, db_module)\n", + "dconstants = getattr(Constants, dclass)(det.detector_type)\n", + "\n", + "print('CalDB Interface: {}'.format(cal_db_interface))\n", + "print('Start time at: ', start)\n", + "print('End time at: ', end)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "parameter_list = combine_lists(*parameters, names = parameter_names)\n", + "print(parameter_list)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "# Retrieve list of meta-data\n", + "constant_versions = []\n", + "constant_parameters = []\n", + "constantBP_versions = []\n", + "\n", + "# Loop over constants\n", + "for c, const in enumerate(constants):\n", + " \n", + " if use_existing != \"\":\n", + " break\n", + " \n", + " # Loop over parameters\n", + " for pars in parameter_list:\n", + " \n", + " if (const in [\"Offset\", \"Noise\", \"SlopesPC\"] or \"DARK\" in const.upper()):\n", + " dcond = Conditions.Dark\n", + " mcond = getattr(dcond, dclass)(**pars)\n", + " else:\n", + " dcond = Conditions.Illuminated\n", + " mcond = getattr(dcond, dclass)(**pars,\n", + " photon_energy=photon_energy)\n", + "\n", + " \n", + " \n", + " print('Request: ', const, 'with paramters:', pars)\n", + " # Request Constant versions for given parameters and module\n", + " data = get_from_db(det,\n", + " getattr(dconstants,\n", + " const)(),\n", + " copy.deepcopy(mcond), None,\n", + " cal_db_interface,\n", + " creation_time=start,\n", + " verbosity=0,\n", + " timeout=cal_db_timeout,\n", + " meta_only=True,\n", + " version_info=True)\n", + " \n", + " if not isinstance(data, list):\n", + " continue\n", + " \n", + " data = sorted(data, key=itemgetter('begin_at'), reverse=True)\n", + " print('Number of retrieved constants: {}'.format(len(data)) )\n", + " \n", + " if const in constantsDark:\n", + " # Request BP constant versions\n", + " dataBP = get_from_db(det,\n", + " getattr(dconstants, \n", + " constantsDark[const])(),\n", + " copy.deepcopy(mcond), None,\n", + " cal_db_interface,\n", + " creation_time=start,\n", + " verbosity=0,\n", + " timeout=cal_db_timeout,\n", + " meta_only=True,\n", + " version_info=True)\n", + " \n", + " if not isinstance(data, list) or not isinstance(dataBP, list):\n", + " continue\n", + " print('Number of retrieved darks: {}'.format(len(dataBP)) )\n", + " found_BPmatch = False\n", + " for d in data:\n", + " # Match proper BP constant version\n", + " # and get constant version within\n", + " # requested time range\n", + " if d is None:\n", + " print('Time or data is not found!')\n", + " continue\n", + "\n", + " dt = dateutil.parser.parse(d['begin_at'])\n", + "\n", + " if (dt.replace(tzinfo=None) > end or \n", + " (nconstants==0 and dt.replace(tzinfo=None) < start)):\n", + " continue\n", + " \n", + " if nconstants>0 and constant_parameters.count(pars)>nconstants-1:\n", + " break\n", + "\n", + " closest_BP = None\n", + " closest_BPtime = None\n", + "\n", + " for dBP in dataBP:\n", + " if dBP is None:\n", + " print(\"Bad pixels are not found!\")\n", + " continue\n", + "\n", + " dt = dateutil.parser.parse(d['begin_at'])\n", + " dBPt = dateutil.parser.parse(dBP['begin_at'])\n", + "\n", + " if dt == dBPt:\n", + " found_BPmatch = True\n", + " else:\n", + "\n", + " if np.abs(dBPt-dt).seconds < (max_time*60):\n", + " if closest_BP is None:\n", + " closest_BP = dBP\n", + " closest_BPtime = dBPt\n", + " else:\n", + " if np.abs(dBPt-dt) < np.abs(closest_BPtime-dt):\n", + " closest_BP = dBP\n", + " closest_BPtime = dBPt\n", + "\n", + " if dataBP.index(dBP) == len(dataBP)-1:\n", + " if closest_BP:\n", + " dBP = closest_BP\n", + " dBPt = closest_BPtime\n", + " found_BPmatch = True\n", + " else:\n", + " print('Bad pixels are not found!')\n", + "\n", + " if found_BPmatch:\n", + " print(\"Found constant {}: begin at {}\".format(const, dt))\n", + " print(\"Found bad pixels at {}\".format(dBPt))\n", + " constantBP_versions.append(dBP)\n", + " constant_versions.append(d)\n", + " constant_parameters.append(copy.deepcopy(pars))\n", + " found_BPmatch = False\n", + " break\n", + " else:\n", + " constant_versions += data\n", + " constant_parameters += [copy.deepcopy(pars)]*len(data)\n", + "\n", + "# Remove dublications\n", + "constant_versions_tmp = []\n", + "constant_parameters_tmp = []\n", + "constantBP_versions_tmp = []\n", + "for i, x in enumerate(constant_versions):\n", + " if x not in constant_versions_tmp:\n", + " constant_versions_tmp.append(x)\n", + " constant_parameters_tmp.append(constant_parameters[i])\n", + " if i<len(constantBP_versions)-1:\n", + " constantBP_versions_tmp.append(constantBP_versions[i])\n", + "constant_versions=constant_versions_tmp\n", + "constantBP_versions=constantBP_versions_tmp\n", + "constant_parameters=constant_parameters_tmp\n", + "\n", + "print('Number of stored constant versions is {}'.format(len(constant_versions)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def get_rebined(a, rebin):\n", + " return a[:,:,0].reshape(\n", + " int(a.shape[0] / rebin[0]),\n", + " rebin[0],\n", + " int(a.shape[1] / rebin[1]),\n", + " rebin[1])\n", + " \n", + "def modify_const(const, data, isBP = False):\n", + " return data\n", + "\n", + "ret_constants = {}\n", + "constant_data = ConstantMetaData()\n", + "constant_BP = ConstantMetaData()\n", + "for i, constant_version in enumerate(constant_versions):\n", + "\n", + " const = constant_version['data_set_name'].split('/')[-2]\n", + " qm = db_module\n", + " \n", + " print(\"constant: {}, module {}\".format(const,qm))\n", + " constant_data.retrieve_from_version_info(constant_version)\n", + " \n", + " for key in separate_plot:\n", + " const = '{}_{}'.format(const, constant_parameters[i][key])\n", + " \n", + " if not const in ret_constants:\n", + " ret_constants[const] = {}\n", + " if not qm in ret_constants[const]:\n", + " ret_constants[const][qm] = []\n", + " \n", + " cdata = constant_data.calibration_constant.data\n", + " ctime = constant_data.calibration_constant_version.begin_at\n", + " \n", + " cdata = modify_const(const, cdata)\n", + " \n", + " if len(constantBP_versions)>0:\n", + " constant_BP.retrieve_from_version_info(constantBP_versions[i])\n", + " cdataBP = constant_BP.calibration_constant.data\n", + " cdataBP = modify_const(const, cdataBP, True)\n", + " \n", + " if cdataBP.shape != cdata.shape:\n", + " print('Wrong bad pixel shape! {}, expected {}'.format(cdataBP.shape, cdata.shape))\n", + " continue\n", + " \n", + " # Apply bad pixel mask\n", + " cdataABP = np.copy(cdata)\n", + " cdataABP[cdataBP > 0] = np.nan\n", + " \n", + " # Create superpixels for constants with BP applied\n", + " cdataABP = get_rebined(cdataABP, spShape)\n", + " toStoreBP = np.nanmean(cdataABP, axis=(1, 3))\n", + " toStoreBPStd = np.nanstd(cdataABP, axis=(1, 3))\n", + "\n", + " # Prepare number of bad pixels per superpixels\n", + " cdataBP = get_rebined(cdataBP, spShape)\n", + " cdataNBP = np.nansum(cdataBP > 0, axis=(1, 3))\n", + " else:\n", + " toStoreBP = 0\n", + " toStoreBPStd = 0\n", + " cdataNBP = 0\n", + "\n", + " # Create superpixels for constants without BP applied\n", + " cdata = get_rebined(cdata, spShape)\n", + " toStoreStd = np.nanstd(cdata, axis=(1, 3))\n", + " toStore = np.nanmean(cdata, axis=(1, 3))\n", + " \n", + " # Convert parameters to dict\n", + " dpar = {p.name: p.value for p in constant_data.detector_condition.parameters}\n", + " \n", + " print(\"Store values in dict\", const, qm, ctime)\n", + " ret_constants[const][qm].append({'ctime': ctime,\n", + " 'nBP': cdataNBP,\n", + " 'dataBP': toStoreBP,\n", + " 'dataBPStd': toStoreBPStd,\n", + " 'data': toStore,\n", + " 'dataStd': toStoreStd,\n", + " 'mdata': dpar}) \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "if use_existing == \"\":\n", + " print('Save data to {}/CalDBAna_{}_{}.h5'.format(out_folder, dclass, db_module))\n", + " save_dict_to_hdf5(ret_constants,\n", + " '{}/CalDBAna_{}_{}.h5'.format(out_folder, dclass, db_module))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if use_existing == \"\":\n", + " fpath = '{}/CalDBAna_{}_*.h5'.format(out_folder, dclass)\n", + "else:\n", + " fpath = '{}/CalDBAna_{}_*.h5'.format(use_existing, dclass)\n", + "\n", + "print('Load data from {}'.format(fpath))\n", + "ret_constants = load_data_from_hdf5(fpath)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Parameters for plotting\n", + "\n", + "keys = {\n", + " 'Mean': ['data', '', 'Mean over pixels'],\n", + " 'std': ['dataStd', '', '$\\sigma$ over pixels'],\n", + " 'MeanBP': ['dataBP', 'Good pixels only', 'Mean over pixels'],\n", + " 'NBP': ['nBP', 'Fraction of BP', 'Fraction of BP'],\n", + " 'stdBP': ['dataBPStd', 'Good pixels only', '$\\sigma$ over pixels'],\n", + " 'stdASIC': ['', '', '$\\sigma$ over ASICs'],\n", + " 'stdCell': ['', '', '$\\sigma$ over Cells'],\n", + "}\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print('Plot calibration constants')\n", + "\n", + "# loop over constat type\n", + "for const, modules in ret_constants.items():\n", + "\n", + " const = const.split(\"_\")\n", + " print('Const: {}'.format(const))\n", + "\n", + " # Loop over modules\n", + " for mod, data in modules.items():\n", + " print(mod)\n", + "\n", + " ctimes = np.array(data[\"ctime\"])\n", + " ctimes_ticks = [x.strftime('%y-%m-%d') for x in ctimes]\n", + "\n", + " if (\"mdata\" in data):\n", + " cmdata = np.array(data[\"mdata\"])\n", + " for i, tick in enumerate(ctimes_ticks):\n", + " ctimes_ticks[i] = ctimes_ticks[i] + \\\n", + " ', V={:1.0f}'.format(cmdata[i]['Sensor Temperature']) + \\\n", + " ', T={:1.0f}'.format(\n", + " cmdata[i]['Integration Time'])\n", + "\n", + " sort_ind = np.argsort(ctimes_ticks)\n", + " ctimes_ticks = list(np.array(ctimes_ticks)[sort_ind])\n", + "\n", + " # Create sorted by data dataset\n", + " rdata = {}\n", + " for key, item in keys.items():\n", + " if item[0] in data:\n", + " rdata[key] = np.array(data[item[0]])[sort_ind]\n", + "\n", + " nTimes = rdata['Mean'].shape[0]\n", + " nPixels = rdata['Mean'].shape[1] * rdata['Mean'].shape[2]\n", + " nBins = nPixels\n", + "\n", + " # Avoid too low values\n", + " if const[0] in [\"Noise\", \"Offset\"]:\n", + " rdata['Mean'][rdata['Mean'] < 0.1] = np.nan\n", + " if 'MeanBP' in rdata:\n", + " rdata['MeanBP'][rdata['MeanBP'] < 0.1] = np.nan\n", + " \n", + " if 'NBP' in rdata:\n", + " rdata['NBP'] = rdata['NBP'].astype(float)\n", + " rdata[\"NBP\"][rdata[\"NBP\"] == (spShape[0] * spShape[1])] = np.nan\n", + " rdata[\"NBP\"] = rdata[\"NBP\"] / spShape[0] / spShape[1] * 100\n", + "\n", + " # Reshape: ASICs over cells for plotting\n", + " pdata = {}\n", + " for key in rdata:\n", + " if len(rdata[key].shape)<3:\n", + " continue\n", + " pdata[key] = rdata[key][:, :, :].reshape(nTimes, nBins).swapaxes(0, 1)\n", + "\n", + " # Plotting\n", + " for key in pdata:\n", + " if len(pdata[key].shape)<2:\n", + " continue\n", + "\n", + " if key == 'NBP':\n", + " unit = '[%]'\n", + " else:\n", + " unit = '[ADU]'\n", + "\n", + " title = '{}, module {}, {}'.format(\n", + " const[0], mod, keys[key][1])\n", + " cb_label = '{}, {} {}'.format(const[0], keys[key][2], unit)\n", + "\n", + " fname = '{}/{}_{}'.format(out_folder, const[0], mod.replace('_', ''))\n", + " for item in const[1:]:\n", + " fname = '{}_{}'.format(fname, item)\n", + " fname = '{}_ASIC_{}.png'.format(fname, key)\n", + " \n", + " vmin,vmax = get_range(pdata[key][::-1].flatten(), plot_range)\n", + " hm_combine(pdata[key][::-1], htype=HMType.mro,\n", + " x_label='Creation Time', y_label='ASIC ID',\n", + " x_ticklabels=ctimes_ticks,\n", + " x_ticks=np.arange(len(ctimes_ticks))+0.3,\n", + " title=title, cb_label=cb_label,\n", + " vmin=vmin, vmax=vmax,\n", + " fname=fname,\n", + " pad=[0.125, 0.125, 0.12, 0.185])\n" + ] + } + ], + "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.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} -- GitLab