From f7111657122ad56d0391f3c04cedb6fd8bf308ac Mon Sep 17 00:00:00 2001
From: Kiana Setoodehnia <kiana.setoodehnia@xfel.eu>
Date: Wed, 29 Apr 2020 14:57:19 +0200
Subject: [PATCH] I was thinking to reduce the bin numbers but I will run into
other issues, so I put it at the original values.
---
.../pnCCD/Characterize_pnCCD_Dark_NBC.ipynb | 30 +-
notebooks/pnCCD/Characterize_pnCCD_Gain.ipynb | 873 +++++----
notebooks/pnCCD/Correct_pnCCD_NBC.ipynb | 1597 ++++++++++-------
xfel_calibrate/notebooks.py | 11 +-
4 files changed, 1567 insertions(+), 944 deletions(-)
diff --git a/notebooks/pnCCD/Characterize_pnCCD_Dark_NBC.ipynb b/notebooks/pnCCD/Characterize_pnCCD_Dark_NBC.ipynb
index 2ebcc78c2..2d96e9ac5 100644
--- a/notebooks/pnCCD/Characterize_pnCCD_Dark_NBC.ipynb
+++ b/notebooks/pnCCD/Characterize_pnCCD_Dark_NBC.ipynb
@@ -48,6 +48,7 @@
"cal_db_timeout = 300000 # timeout on caldb requests\n",
"db_output = False # if True, the notebook sends dark constants to the calibration database\n",
"local_output = True # if True, the notebook saves dark constants locally\n",
+ "creation_time = \"\" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC.00 e.g. 2019-07-04 11:02:41.00\n",
"\n",
"number_dark_frames = 0 # number of images to be used, if set to 0 all available images are used\n",
"chunkSize = 100 # number of images to read per chunk\n",
@@ -79,6 +80,7 @@
"source": [
"import os\n",
"import copy\n",
+ "import datetime\n",
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"\n",
@@ -169,18 +171,29 @@
"filename = fp_path.format(sequence)\n",
"h5path = h5path.format(karabo_id, receiver_id)\n",
"\n",
- "creation_time = None\n",
- "if use_dir_creation_date:\n",
+ "# Run's creation time:\n",
+ "if creation_time:\n",
+ " try:\n",
+ " creation_time = datetime.datetime.strptime(creation_time, '%Y-%m-%d %H:%M:%S.%f')\n",
+ " except Exception as e:\n",
+ " print(f\"creation_time value error: {e}.\" \n",
+ " \"Use same format as YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00/n\")\n",
+ " creation_time = None\n",
+ " print(\"Given creation time wont be used.\")\n",
+ "else:\n",
+ " creation_time = None\n",
+ "\n",
+ "if not creation_time and use_dir_creation_date:\n",
" creation_time = get_dir_creation_date(in_folder, run)\n",
+ "\n",
+ "print(f\"Creation time: {creation_time}\")\n",
" \n",
"cal_db_interface = get_random_db_interface(cal_db_interface)\n",
"print('Calibration database interface: {}'.format(cal_db_interface))\n",
"print(\"Sending constants to the calibration database: {}\".format(db_output))\n",
"print(\"HDF5 path to data: {}\".format(h5path))\n",
"print(\"Reading data from: {}\".format(filename))\n",
- "print(\"Run number: {}\".format(run))\n",
- "if creation_time:\n",
- " print(\"Creation time: {}\".format(creation_time))"
+ "print(\"Run number: {}\".format(run))"
]
},
{
@@ -1019,6 +1032,13 @@
" f\"• gain_setting: {gain}\\n• temperature: {temperature_k}\\n\"\n",
" f\"• creation_time: {creation_time}\\n\")"
]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
}
],
"metadata": {
diff --git a/notebooks/pnCCD/Characterize_pnCCD_Gain.ipynb b/notebooks/pnCCD/Characterize_pnCCD_Gain.ipynb
index 8ed501a08..6cb9c9ef1 100644
--- a/notebooks/pnCCD/Characterize_pnCCD_Gain.ipynb
+++ b/notebooks/pnCCD/Characterize_pnCCD_Gain.ipynb
@@ -4,11 +4,15 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# pnCCD Gain Characterization ##\n",
+ "# pnCCD Gain Characterization #\n",
"\n",
- "Authors: DET Group, Version 1.0\n",
+ "Authors: DET Group, modified by Kiana Setoodehnia on March 2020 - Version 2.0\n",
"\n",
- "The following notebook provides gain characterization for the pnCCD. It relies on data previously having been corrected using the MDC interface"
+ "The following notebook provides gain characterization for the pnCCD. It relies on data which are previously corrected using the Meta Data Catalog web service interface or by running the Correct_pnCCD_NBC.ipynb notebook. The corrections which are applied by the web service or the aforementioned notebook are as follows:\n",
+ "\n",
+ "- offset correction\n",
+ "- common mode correction\n",
+ "- split pattern classification"
]
},
{
@@ -18,56 +22,64 @@
"ExecuteTime": {
"end_time": "2018-12-06T15:54:23.218849Z",
"start_time": "2018-12-06T15:54:23.166497Z"
- },
- "collapsed": true
+ }
},
"outputs": [],
"source": [
+ "cluster_profile = \"noDB\" # ipcluster profile to use\n",
+ "in_folder = \"/gpfs/exfel/exp/SQS/201930/p900075/proc\" # input folder\n",
+ "out_folder = '/gpfs/exfel/exp/SQS/201930/p900075/proc' # output folder\n",
+ "run = 365 # which run to read data from\n",
+ "sequences = [-1] # sequences to correct, set to -1 for all, range allowed\n",
"\n",
- "in_folder = \"/gpfs/exfel/exp/SQS/201921/p002430/proc/\" # input folder, required\n",
- "out_folder = '/gpfs/exfel/data/scratch/xcal/test/' # output folder, required\n",
+ "karabo_da = 'PNCCD01' # data aggregators\n",
+ "karabo_id = \"SQS_NQS_PNCCD1MP\" # karabo prefix of PNCCD devices\n",
+ "receiver_id = \"PNCCD_FMT-0\" # inset for receiver devices\n",
"path_template = 'CORR-R{:04d}-PNCCD01-S{{:05d}}.h5' # the template to use to access data\n",
"path_template_seqs = \"{}/r{:04d}/*PNCCD01-S*.h5\"\n",
- "run = 747 # which run to read data from, required\n",
- "cluster_profile = \"noDB\" # ipcluster profile to use\n",
- "sigma_noise = 10. # Pixel exceeding 'sigmaNoise' * noise value in that pixel will be masked\n",
- "h5path = '/INSTRUMENT/SQS_NQS_PNCCD1MP/CAL/PNCCD_FMT-0:output/data/' # path in the HDF5 file the data is at\n",
- "multi_iteration = False # use multiple iterations\n",
- "use_dir_creation_date = True\n",
+ "h5path = '/INSTRUMENT/{}/CAL/{}:output/data/' # path to data in the HDF5 file \n",
+ "\n",
+ "cpuCores = 8\n",
+ "use_dir_creation_date = True # this is needed to obtain creation time of the run\n",
+ "overwrite = True # keep this as True to not overwrite the output\n",
+ "sequences_per_node = 1\n",
+ "chunkSize = 100 # number of images to read per chunk\n",
+ "run_parallel = True\n",
+ "db_output = False # if True, the notebook injects dark constants into the calibration database\n",
+ "local_output = True # if True, the notebook saves dark constants locally\n",
+ "\n",
+ "cal_db_interface = \"tcp://max-exfl016:8015\" # calibration DB interface to use\n",
+ "cal_db_timeout = 300000 # timeout on caldb requests\n",
+ "creation_time = \"\" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC.00 e.g. 2019-07-04 11:02:41.00\n",
+ "\n",
+ "# pnCCD parameters:\n",
"fix_temperature = 233.\n",
- "gain = 0\n",
+ "gain = 1\n",
"bias_voltage = 300\n",
"integration_time = 70\n",
- "cal_db_interface = \"tcp://max-exfl016:8015\" # 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\n",
- "sequences_per_node = 1\n",
- "limit_images = 0\n",
- "time_offset_days = 0\n",
- "cpuCores = 8\n",
- "\n",
+ "photon_energy = 1.6 # Al fluorescence in keV\n",
"\n",
- "def balance_sequences(in_folder, run, sequences, sequences_per_node, path_template_seqs):\n",
- " import glob\n",
- " import re\n",
- " import numpy as np\n",
- " if sequences[0] == -1:\n",
- " sequence_files = glob.glob(path_template_seqs.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",
- " 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"
+ "def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da):\n",
+ " from xfel_calibrate.calibrate import balance_sequences as bs\n",
+ " return bs(in_folder, run, sequences, sequences_per_node, karabo_da)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# On the singles spectrum (uploaded in the middle of this notebook), the ADU values correspoding to the boundaries\n",
+ "# of the first peak region are used as cti_limit_low and cti_limit_high:\n",
+ "if gain == 1:\n",
+ " cti_limit_low = 3000 # lower limit of cti\n",
+ " cti_limit_high = 10000 # higher limit of cti\n",
+ " max_points = 100000 # maximum data value\n",
+ "elif gain == 64:\n",
+ " cti_limit_low = 50 \n",
+ " cti_limit_high = 170 \n",
+ " max_points = 500000 "
]
},
{
@@ -77,11 +89,28 @@
"ExecuteTime": {
"end_time": "2018-12-06T15:54:23.455376Z",
"start_time": "2018-12-06T15:54:23.413579Z"
- },
- "collapsed": true
+ }
},
"outputs": [],
"source": [
+ "import numpy as np\n",
+ "import h5py\n",
+ "import matplotlib.pyplot as plt\n",
+ "import iminuit as im\n",
+ "from iminuit import Minuit\n",
+ "from IPython.display import display, Markdown\n",
+ "import copy\n",
+ "import glob\n",
+ "import os\n",
+ "from prettytable import PrettyTable\n",
+ "import datetime\n",
+ "from datetime import timedelta\n",
+ "from mpl_toolkits.axes_grid1 import ImageGrid, AxesGrid\n",
+ "from functools import partial\n",
+ "import matplotlib\n",
+ "import warnings\n",
+ "warnings.filterwarnings('ignore')\n",
+ "\n",
"import XFELDetAna.xfelprofiler as xprof\n",
"\n",
"profiler = xprof.Profiler()\n",
@@ -89,41 +118,20 @@
"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.detectors.fastccd import readerh5 as fastccdreaderh5\n",
"from XFELDetAna.plotting.util import prettyPlotting\n",
"prettyPlotting=True\n",
"from XFELDetAna.xfelreaders import ChunkReader\n",
- "\n",
- "import numpy as np\n",
- "import h5py\n",
- "import matplotlib.pyplot as plt\n",
- "from iminuit import Minuit\n",
- "import copy\n",
- "import os\n",
- "\n",
- "from prettytable import PrettyTable\n",
- "\n",
"from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions\n",
- "from cal_tools.tools import get_dir_creation_date\n",
- "\n",
- "import datetime\n",
- "from datetime import timedelta\n",
+ "from iCalibrationDB.detectors import DetectorTypes\n",
+ "from cal_tools.tools import get_dir_creation_date, save_const_to_h5\n",
"\n",
"%matplotlib inline\n",
"\n",
"if sequences[0] == -1:\n",
- " sequences = None\n",
- "\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)))])"
+ " sequences = None"
]
},
{
@@ -133,33 +141,65 @@
"ExecuteTime": {
"end_time": "2018-12-06T15:54:23.679069Z",
"start_time": "2018-12-06T15:54:23.662821Z"
- },
- "collapsed": false
+ }
},
"outputs": [],
"source": [
- "\n",
- "x = 1024 # rows of the FastCCD to analyze in FS mode \n",
- "y = 1024 # columns of the FastCCD to analyze in FS mode \n",
- "sensorSize = [x, y]\n",
- " \n",
+ "# Calibration Database Settings, and Some Initial Run Parameters & Paths:\n",
+ "display(Markdown('### Initial Settings and Paths'))\n",
+ "pixels_x = 1024 # rows of pnCCD in pixels \n",
+ "pixels_y = 1024 # columns of pnCCD in pixels\n",
+ "print(f\"pnCCD size is: {pixels_x}x{pixels_y} pixels.\")\n",
+ "print(f'Calibration database interface to be used: {cal_db_interface}')\n",
+ "\n",
+ "proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]\n",
+ "file_loc =f'Proposal: {proposal}, Run: {run}'\n",
+ "print(file_loc)\n",
+ "\n",
+ "# Paths to the data:\n",
"ped_dir = \"{}/r{:04d}\".format(in_folder, run)\n",
- "fp_name = path_template.format(run)\n",
- "\n",
+ "fp_name = path_template.format(run, karabo_da)\n",
+ "fp_path = '{}/{}'.format(ped_dir, fp_name)\n",
+ "h5path = h5path.format(karabo_id, receiver_id)\n",
+ "print(\"HDF5 path to data: {}\\n\".format(h5path))\n",
"\n",
- "creation_time = None\n",
- "if use_dir_creation_date:\n",
+ "# Run's creation time:\n",
+ "if creation_time:\n",
+ " try:\n",
+ " creation_time = datetime.datetime.strptime(creation_time, '%Y-%m-%d %H:%M:%S.%f')\n",
+ " except Exception as e:\n",
+ " print(f\"creation_time value error: {e}.\" \n",
+ " \"Use same format as YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00/n\")\n",
+ " creation_time = None\n",
+ " print(\"Given creation time wont be used.\")\n",
+ "else:\n",
+ " creation_time = None\n",
+ "\n",
+ "if not creation_time and use_dir_creation_date:\n",
" creation_time = get_dir_creation_date(in_folder, run)\n",
"\n",
- "fp_path = '{}/{}'.format(ped_dir, fp_name)\n",
+ "print(f\"Creation time: {creation_time}\")\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()))\n",
- " \n",
- "out_folder = \"{}/r{:04d}\".format(out_folder, run)"
+ "# Reading all sequences of the run:\n",
+ "file_list = []\n",
+ "total_sequences = 0\n",
+ "fsequences = []\n",
+ "\n",
+ "if sequences is None:\n",
+ " file_list = glob.glob(fp_path.format(0).replace('00000', '*'))\n",
+ " file_list = sorted(file_list, key = lambda x: (len (x), x))\n",
+ " total_sequences = len(file_list)\n",
+ " fsequences = range(total_sequences)\n",
+ "else:\n",
+ " for seq in sequences:\n",
+ " abs_entry = fp_path.format(seq)\n",
+ " if os.path.isfile(abs_entry):\n",
+ " file_list.append(abs_entry)\n",
+ " total_sequences += 1\n",
+ " fsequences.append(seq)\n",
+ "\n",
+ "sequences = fsequences \n",
+ "print(f\"This run has a total number of {total_sequences} sequences.\\n\")"
]
},
{
@@ -169,29 +209,14 @@
"ExecuteTime": {
"end_time": "2018-12-06T15:54:23.913269Z",
"start_time": "2018-12-06T15:54:23.868910Z"
- },
- "collapsed": true
+ }
},
"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 = [512, 512]\n",
- "commonModeAxisR = 'row'#Axis along which common mode will be calculated\n",
- "run_parallel = True\n",
- "profile = False\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"
+ "display(Markdown('### List of Files to be Processed'))\n",
+ "print(\"Reading data from the following files:\")\n",
+ "for i in range(len(file_list)):\n",
+ " print(file_list[i])"
]
},
{
@@ -201,68 +226,29 @@
"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",
+ "# Sensor size and block size definitions (important for common mode and other corrections):\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).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).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": "markdown",
- "metadata": {},
- "source": [
- "### List of files to be processed ###"
+ "sensorSize = [pixels_x, pixels_y]\n",
+ "blockSize = [sensorSize[0]//2, sensorSize[1]//2] # Sensor area will be analysed according to blocksize\n",
+ "xcal.defaultBlockSize = blockSize\n",
+ "memoryCells = 1 # pnCCD has 1 memory cell"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2018-12-06T18:43:39.776018Z",
- "start_time": "2018-12-06T18:43:39.759185Z"
- },
- "collapsed": false
- },
+ "metadata": {},
"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\"]))) "
+ "# Output Folder Creation:\n",
+ "out_folder = \"{}/r{:04d}\".format(out_folder, run)\n",
+ "os.makedirs(out_folder, exist_ok=True)\n",
+ "if not overwrite:\n",
+ " raise AttributeError(\"Output path exists! Exiting\") "
]
},
{
@@ -271,21 +257,17 @@
"source": [
"### Reading in already corrected data ###\n",
"\n",
- "In the following we read in already corrected data and separte double event back out\n",
- "such as to track split event ratios.\n",
- "\n",
- "Additionally, data structures for CTI and relative gain determination are created"
+ "In the following, we read in already corrected data and separate double events to track the ratios of split events. Additionally, data structures for CTI and relative gain determination are created."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": true
- },
+ "metadata": {},
"outputs": [],
"source": [
- "## pattern kernels as defined in pyDetLib\n",
+ "# pattern kernels as defined in pyDetLib:\n",
+ "# possibly at pydetlib/lib/src/XFELDetAna/util/pattern_kernels.py\n",
"\n",
"kernel = {}\n",
"\n",
@@ -336,8 +318,6 @@
" idx = np.roll(pdat, 1, axis = 1) == p\n",
" p2 = alldat[idx]\n",
" \n",
- " \n",
- " \n",
" r = np.zeros((p1.shape[0], 2), np.float)\n",
" r[:,0] = np.abs(p1)-np.abs(p2)\n",
" r[:,1] = np.abs(p2)\n",
@@ -352,12 +332,11 @@
"ExecuteTime": {
"end_time": "2018-12-06T16:10:55.917179Z",
"start_time": "2018-12-06T16:09:01.603633Z"
- },
- "collapsed": true
+ }
},
"outputs": [],
"source": [
- "# actual data reading happens here\n",
+ "# actual data reading happens here:\n",
"\n",
"# separated doubles will be held in these lists\n",
"up_doubles = [] # up major partner\n",
@@ -365,37 +344,38 @@
"left_doubles = [] # left major partner\n",
"right_doubles = [] # right major partner\n",
"\n",
- "x = 1024\n",
- "y = 1024\n",
+ "# A struct in the C programming language is a composite data type declaration that defines a physically grouped \n",
+ "# list of variables under one name in a block of memory, allowing the different variables to be accessed via a \n",
+ "# single pointer or by the struct declared name which returns the same address. \n",
"\n",
- "# structs for CTI and gain evaluation\n",
- "# we work by quadrant\n",
- "dstats = {'upper left': {'coords': [0, x//2, 0, y//2], # coordinates of this quadrant\n",
+ "# structs for CTI and gain evaluation:\n",
+ "# we work by quadrant:\n",
+ "dstats = {'upper left': {'coords': [0, pixels_x//2, 0, pixels_y//2], # coordinates of this quadrant\n",
" 'row_coords': [], # will hold our data's row coordinates\n",
" 'col_coords': [], # will hold our data's column coordinates\n",
" 'col_dir': 1, # 1 if columns read out to left, -1 if read out to right\n",
" 'vals': [], # will hold our data values\n",
" },\n",
- " 'lower left': {'coords':[0, x//2, y//2, y],\n",
+ " 'lower left': {'coords':[0, pixels_x//2, pixels_y//2, pixels_y],\n",
" 'row_coords': [],\n",
" 'col_coords': [],\n",
" 'col_dir': 1,\n",
" 'vals': [],\n",
" },\n",
- " 'upper right': {'coords':[x//2, x, 0, y//2],\n",
+ " 'upper right': {'coords':[pixels_x//2, pixels_x, 0, pixels_y//2],\n",
" 'row_coords': [],\n",
" 'col_coords': [],\n",
" 'col_dir': -1,\n",
" 'vals': [],\n",
" },\n",
- " 'lower right': {'coords': [x//2, x, y//2, y],\n",
+ " 'lower right': {'coords': [pixels_x//2, pixels_x, pixels_y//2, pixels_y],\n",
" 'row_coords': [],\n",
" 'col_coords': [],\n",
" 'col_dir': -1,\n",
" 'vals': [],\n",
" },\n",
" \n",
- "}\n",
+ "} \n",
"\n",
"# structure that identifies and holds pattern statistics\n",
"pstats = {'singles': {'p': 100,},\n",
@@ -424,24 +404,24 @@
" patterns = infile[h5path+\"/patterns\"][()]\n",
" bpix = infile[h5path+\"/mask\"][()]\n",
" \n",
- " # separate out doubles\n",
+ " # separating out doubles\n",
" up_doubles.append(separateDouble(patterns, data, 202))\n",
" down_doubles.append(separateDouble(patterns, data, 203))\n",
" right_doubles.append(separateDouble(patterns, data, 201))\n",
" left_doubles.append(separateDouble(patterns, data, 200))\n",
" \n",
- " # create pattern statistics\n",
+ " # creating pattern statistics\n",
" for pstr, entry in pstats.items():\n",
" cpat = patterns == entry[\"p\"]\n",
" entry[\"counts\"] = np.count_nonzero(cpat)\n",
" hist = entry.get(\"hist\", 0)\n",
" hist += np.sum(cpat, axis=0)\n",
- " entry[\"hist\"] = hist\n",
+ " entry[\"hist\"] = hist \n",
" \n",
- " # for relative norms\n",
+ " # for normalization of relative occurences\n",
" tpats += np.count_nonzero((patterns > 0) & (patterns < 1000))\n",
" \n",
- " # create coordinates needed for CTI and relative gain fitting\n",
+ " # creating coordinates needed for CTI and relative gain fitting\n",
" xv, yv = np.meshgrid(np.arange(data.shape[0]), np.arange(data.shape[1]))\n",
" xv = np.repeat(xv[...,None], data.shape[2], axis=2)\n",
" yv = np.repeat(yv[...,None], data.shape[2], axis=2)\n",
@@ -449,17 +429,20 @@
" # we take all single values\n",
" all_singles = np.zeros_like(data)\n",
" all_singles[(patterns == 100) | (patterns == 101)] = data[(patterns == 100) | (patterns == 101)]\n",
+ " #all_singles[bpix != 0] = np.nan\n",
" \n",
" # and then save data and coordinates quadrant-wise\n",
" for key, entry in dstats.items():\n",
" co = entry[\"coords\"]\n",
" cd = entry[\"col_dir\"]\n",
" asd = all_singles[:, co[2]:co[3], co[0]:co[1]]\n",
- " # take column direction into account - right hemisphere reads to right, left to left\n",
+ " # taking column direction into account - right hemisphere reads to right, left hemisphere reads to left\n",
" if cd == 1:\n",
- " yv, xv = np.meshgrid(np.arange(asd.shape[1]), np.arange(asd.shape[2])+co[2]) # row counting downwards\n",
+ " yv, xv = np.meshgrid(np.arange(asd.shape[1]), np.arange(asd.shape[2])+co[2]) # row counting \n",
+ " # downwards\n",
" else:\n",
- " yv, xv = np.meshgrid(np.arange(asd.shape[1], 0, -1), np.arange(asd.shape[2])+co[2]) # row counting downwards\n",
+ " yv, xv = np.meshgrid(np.arange(asd.shape[1], 0, -1), np.arange(asd.shape[2])+co[2]) # row counting \n",
+ " # downwards\n",
" xv = np.repeat(xv[None,...], data.shape[0], axis=0)\n",
" yv = np.repeat(yv[None,...], data.shape[0], axis=0)\n",
" entry['row_coords'].append(xv[asd > 0].flatten())\n",
@@ -479,20 +462,11 @@
" entry['vals'] = np.concatenate(entry['vals'])"
]
},
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": true
- },
- "outputs": [],
- "source": []
- },
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Pattern Statistics ##\n",
+ "### Pattern Statistics ###\n",
"\n",
"Relative occurrences are normalized to non-cluster events"
]
@@ -500,13 +474,11 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": false
- },
+ "metadata": {},
"outputs": [],
"source": [
"table = PrettyTable()\n",
- "table.field_names = [\"Pattern type\", \"Absolute #\", \"Relative Occurence\"]\n",
+ "table.field_names = [\"Pattern type\", \"Absolute Number\", \"Relative Occurence\"]\n",
"for key, entry in pstats.items():\n",
" table.add_row([key, entry[\"counts\"], \"{:0.2f}\".format(float(entry[\"counts\"])/tpats)])\n",
" \n",
@@ -516,12 +488,11 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": false
- },
+ "metadata": {},
"outputs": [],
"source": [
- "from mpl_toolkits.axes_grid1 import ImageGrid\n",
+ "display(Markdown('### Event Patterns'))\n",
+ "\n",
"fig = plt.figure(1, (15., 15.))\n",
"grid = ImageGrid(fig, 111, # similar to subplot(111)\n",
" nrows_ncols=(5, 4), # creates 2x2 grid of axes\n",
@@ -531,63 +502,66 @@
"i = 0\n",
"for key, entry in pstats.items():\n",
" d = entry[\"hist\"]\n",
- " \n",
- " grid[i].imshow(d, vmin=0, vmax=2.*np.mean(d[d != 0]), cmap=\"gray_r\") # The AxesGrid object work as a list of axes.\n",
+ " # The AxesGrid object works as a list of axes.\n",
+ " grid[i].imshow(d, vmin=0, vmax=2.*np.mean(d[d != 0]), cmap=\"gray_r\") \n",
+ " grid[i].set_ylabel(\"Row\")\n",
" grid[i].annotate(key, (20, -30), color=\"red\", annotation_clip=False)\n",
" i += 1\n",
- " if i == 3:\n",
- " i += 1\n"
+ " if i == 3: # For the empty panel\n",
+ " i += 1\n",
+ " if i == 16: # I actually don't know how this is working and why only one panel has x-label\n",
+ " grid[i].set_xlabel(\"Column\")"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": false
- },
+ "metadata": {},
"outputs": [],
"source": [
- "fig = plt.figure(figsize=(8,8))\n",
+ "fig = plt.figure(figsize=(8, 8))\n",
"ax = fig.add_subplot(111)\n",
- "ax.scatter(np.abs(up_doubles[:,0]), np.abs(up_doubles[:,1]), 0.5, alpha=0.1, color='b', label=\"up doubles\")\n",
- "ax.scatter(np.abs(down_doubles[:,1]), np.abs(down_doubles[:,0]), 0.5, alpha=0.1, color='g', label=\"down doubles\")\n",
- "ax.semilogx()\n",
- "ax.semilogy()\n",
- "ax.set_xlim(1000, 50000)\n",
- "ax.set_ylim(1000, 50000)"
+ "ax.scatter(np.abs(up_doubles[:,0]), np.abs(up_doubles[:,1]), 0.5, alpha=0.5, color='b', label=\"up doubles\")\n",
+ "ax.scatter(np.abs(down_doubles[:,1]), np.abs(down_doubles[:,0]), 0.5, alpha=0.5, color='g', label=\"down doubles\")\n",
+ "if gain == 1:\n",
+ " ax.semilogx()\n",
+ " ax.semilogy()\n",
+ " ax.set_xlim(1000, 50000)\n",
+ " ax.set_ylim(1000, 50000)\n",
+ "#ax.set_xlabel('Doubles')\n",
+ "#ax.set_ylabel('Doubles')\n",
+ "ax.legend()\n",
+ "plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": false
- },
+ "metadata": {},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(8,8))\n",
"ax = fig.add_subplot(111)\n",
- "ax.scatter(np.abs(left_doubles[:,0]), np.abs(left_doubles[:,1]), 0.5, alpha=0.1, color='b', label=\"left doubles\")\n",
- "ax.scatter(np.abs(right_doubles[:,1]), np.abs(right_doubles[:,0]), 0.5, alpha=0.1, color='g', label=\"right doubles\")\n",
- "ax.semilogx()\n",
- "ax.semilogy()\n",
- "ax.set_xlim(1000, 50000)\n",
- "ax.set_ylim(1000, 50000)"
+ "ax.scatter(np.abs(left_doubles[:,0]), np.abs(left_doubles[:,1]), 0.5, alpha=0.5, color='b', label=\"left doubles\")\n",
+ "ax.scatter(np.abs(right_doubles[:,1]), np.abs(right_doubles[:,0]), 0.5, alpha=0.5, color='g', label=\"right doubles\")\n",
+ "if gain == 1:\n",
+ " ax.semilogx()\n",
+ " ax.semilogy()\n",
+ " ax.set_xlim(1000, 50000)\n",
+ " ax.set_ylim(1000, 50000) \n",
+ "#ax.set_xlabel('Doubles')\n",
+ "#ax.set_ylabel('Doubles')\n",
+ "ax.legend()\n",
+ "plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": false
- },
+ "metadata": {},
"outputs": [],
"source": [
- "import matplotlib\n",
- "from mpl_toolkits.axes_grid1 import AxesGrid\n",
- "cti_limit_low = 3000\n",
- "cti_limit_high = 10000\n",
- "max_points = 100000\n",
+ "display(Markdown('### Display of Signal along the Columns'))\n",
"\n",
"fig = plt.figure(figsize=(15., 15.))\n",
"grid = AxesGrid(fig, 111, # similar to subplot(111)\n",
@@ -609,7 +583,8 @@
"\n",
" xax = ccoords[~idx]\n",
" yax = rcoords[~idx]\n",
- " grid[i].scatter(xax, avals[~idx], 1, color=\"grey\", alpha=0.2)\n",
+ " grid[i].scatter(xax, avals[~idx], 1, color=\"grey\", alpha=0.2) # These are the data points which ar outside the\n",
+ " # valid CTI limits\n",
"\n",
"\n",
" xax = ccoords[idx]\n",
@@ -623,10 +598,10 @@
" colors = cmap(yax/np.max(yax))\n",
" grid[i].scatter(xax, pvals, 1, color=colors)\n",
" grid[i].scatter(unique_coords, mean_signal, 3, color='tomato')\n",
- " grid[i].set_ylabel(\"Signal value (ADU)\")\n",
- " grid[i].set_xlabel(\"Column coordinate\")\n",
+ " grid[i].set_ylabel(\"Signal Value (ADU)\")\n",
+ " grid[i].set_xlabel(\"Column Coordinate\")\n",
" grid[i].annotate(key, (0.1, 1.02), xycoords=\"axes fraction\", color=\"red\", annotation_clip=False)\n",
- " if entry[\"col_dir\"] == -1: # columns counting from the right\n",
+ " if entry[\"col_dir\"] == -1: # columns reading out to the right\n",
" grid[i].set_xlim(512, 0)\n",
" l = grid[i].set_ylim(0, 1.25*cti_limit_high)\n",
" ai += 1"
@@ -635,12 +610,11 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": false
- },
+ "metadata": {},
"outputs": [],
"source": [
- "\n",
+ "display(Markdown('### Display of Signal along the Rows'))\n",
+ " \n",
"fig = plt.figure(figsize=(15., 15.))\n",
"grid = AxesGrid(fig, 111, # similar to subplot(111)\n",
" nrows_ncols=(2, 2), # creates 2x2 grid of axes\n",
@@ -660,7 +634,8 @@
" \n",
" yax = ccoords[~idx]\n",
" xax = rcoords[~idx]\n",
- " grid[i].scatter(avals[~idx], xax, 1, color=\"grey\", alpha=0.2)\n",
+ " grid[i].scatter(avals[~idx], xax, 1, color=\"grey\", alpha=0.2) # These are the data points which ar outside the\n",
+ " # valid CTI limits\n",
"\n",
" yax = ccoords[idx]\n",
" xax = rcoords[idx]\n",
@@ -673,31 +648,129 @@
" colors = cmap(yax/np.max(yax))\n",
" grid[i].scatter(pvals, xax, 1, color=colors)\n",
" grid[i].scatter(mean_signal, unique_coords, 3, color='tomato')\n",
- " grid[i].set_xlabel(\"Signal value (ADU)\")\n",
- " grid[i].set_ylabel(\"Row coordinate\")\n",
+ " grid[i].set_xlabel(\"Signal Value (ADU)\")\n",
+ " grid[i].set_ylabel(\"Row Coordinate\")\n",
" grid[i].annotate(key, (0.1, 1.02), xycoords=\"axes fraction\", color=\"red\", annotation_clip=False)\n",
" grid[i].set_ylim(co[3], co[2])\n",
" l = grid[i].set_xlim(0, 1.25*cti_limit_high)\n",
" ai += 1"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Fitting is performed in this cell for just a randomly selected sample row as an example:\n",
+ "\n",
+ "def show_fit(row_example):\n",
+ " \n",
+ " if row_example > 512 or row_example < 0:\n",
+ " print(\"Invalid row number. Try again.\")\n",
+ " else:\n",
+ " max_points = -1\n",
+ " limit_cols = [0, 512]\n",
+ "\n",
+ " for key, entry in dstats.items():\n",
+ "\n",
+ " rcoords = entry[\"row_coords\"][:max_points]\n",
+ " ccoords = entry[\"col_coords\"][:max_points]\n",
+ " avals = entry[\"vals\"][:max_points]\n",
+ " co = entry[\"coords\"]\n",
+ " idx = (avals > cti_limit_low) & (avals < cti_limit_high) & (ccoords >= limit_cols[0]) & \\\n",
+ " (ccoords <= limit_cols[1])\n",
+ " xax = rcoords[idx]-co[2]\n",
+ " yax = ccoords[idx] # columns\n",
+ " avalsf = avals[idx]\n",
+ " min_data = 5\n",
+ "\n",
+ " ms = np.zeros(sensorSize[0]//2, np.float)\n",
+ " bs = np.zeros(sensorSize[0]//2, np.float)\n",
+ " mErr = np.zeros(sensorSize[0]//2, np.float)\n",
+ " bErr = np.zeros(sensorSize[0]//2, np.float)\n",
+ "\n",
+ " axes = None\n",
+ " fig = None\n",
+ " fig, ax = plt.subplots(1, sharex=True, figsize=(10, 5))\n",
+ "\n",
+ " def linFunc(x, *p):\n",
+ " return p[1] * x + p[0]\n",
+ "\n",
+ " for i in [row_example]:\n",
+ " idx = (xax == i)\n",
+ " zm = avalsf[idx]\n",
+ " zx = yax[idx]\n",
+ "\n",
+ " def linWrap(m, b):\n",
+ " return np.sum(((linFunc(zx, b, -m) - zm)) ** 2)\n",
+ "\n",
+ " pparm = dict(throw_nan=False, pedantic=False, print_level=0)\n",
+ " pparm[\"m\"] = 0\n",
+ " if gain == 1:\n",
+ " pparm[\"limit_m\"] = (0.9, 1)\n",
+ " elif gain == 64:\n",
+ " pparm[\"limit_m\"] = (0, 1)\n",
+ " pparm[\"b\"] = np.nanmean(zm)\n",
+ "\n",
+ " mini = Minuit(linWrap, **pparm)\n",
+ " fmin = mini.migrad()\n",
+ "\n",
+ " if fmin[0]['is_valid'] and fmin[0]['has_covariance'] and (min_data is None or zx.size > min_data):\n",
+ "\n",
+ " ferr = mini.minos(maxcall = 10000)\n",
+ " res = mini.values\n",
+ "\n",
+ " ms[i] = res[\"m\"]\n",
+ " bs[i] = res[\"b\"]\n",
+ "\n",
+ " mErr[i] = max(-ferr['m']['lower'], ferr['m']['upper'])\n",
+ " bErr[i] = max(-ferr['b']['lower'], ferr['b']['upper'])\n",
+ "\n",
+ " # I am commenting out these two lines because they seem to reduce the number of events in each \n",
+ " # panel:\n",
+ " #zx, uidx = np.unique(zx, return_index=True)\n",
+ " #zm = zm[uidx]\n",
+ "\n",
+ " yt = linFunc(zx, res[\"b\"], res[\"m\"])\n",
+ "\n",
+ " ax.scatter(zx, zm, 3) # point size = 3\n",
+ " ax.plot(zx, -ms[i] * zx + bs[i], color='r',\n",
+ " label = 'CTI fit for row {} in {} quadrant'.format(row_example, key))\n",
+ " ax.annotate('valid: {}'.format(fmin[0]['is_valid']), \n",
+ " xy=(0.2, 0.3), xycoords='axes fraction')\n",
+ " ax.set_xlabel(\"Column Coordinate\")\n",
+ " ax.set_ylabel(\"Singles Events (between CTI_low and CTI_high)\")\n",
+ " ax.set_title(\"Events are only for One Row\")\n",
+ " ax.set_xlim(0, 512)\n",
+ " ax.legend()"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
"metadata": {
- "collapsed": true,
"scrolled": false
},
"outputs": [],
"source": [
- "from functools import partial\n",
- "from iminuit import Minuit\n",
- "\n",
- "\n",
+ "display(Markdown('### Example of CTI Fit for a Specific Row'))\n",
+ "fig = show_fit(365)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "# Fitting is performed in this cell for the entire detector to get the CTI and relative gain:\n",
"\n",
"max_points = -1\n",
"limit_cols = [0, 512]\n",
- "\n",
+ "counter = 0 # Total number of bad columns (over all quadrants) for which the fit does not converge\n",
"\n",
"for key, entry in dstats.items():\n",
" \n",
@@ -705,7 +778,8 @@
" ccoords = entry[\"col_coords\"][:max_points]\n",
" avals = entry[\"vals\"][:max_points]\n",
" co = entry[\"coords\"]\n",
- " idx = (avals > cti_limit_low) & (avals < cti_limit_high) & (ccoords >= limit_cols[0]) & (ccoords <= limit_cols[1])\n",
+ " idx = (avals > cti_limit_low) & (avals < cti_limit_high) & (ccoords >= limit_cols[0]) & \\\n",
+ " (ccoords <= limit_cols[1])\n",
" xax = rcoords[idx]-co[2]\n",
" yax = ccoords[idx]\n",
" avalsf = avals[idx]\n",
@@ -715,37 +789,32 @@
" bs = np.zeros(sensorSize[0]//2, np.float)\n",
" mErr = np.zeros(sensorSize[0]//2, np.float)\n",
" bErr = np.zeros(sensorSize[0]//2, np.float)\n",
- " residuals = np.zeros((sensorSize[0]//2, sensorSize[0]//2), np.float)\n",
- " \n",
- " rtest = np.arange(sensorSize[0]//2)\n",
"\n",
" def linFunc(x, *p):\n",
" return p[1] * x + p[0]\n",
" \n",
- " \n",
- " \n",
" for i in range(sensorSize[0]//2):\n",
" idx = (xax == i)\n",
" zm = avalsf[idx]\n",
" zx = yax[idx]\n",
- "\n",
" \n",
" def linWrap(m, b):\n",
" return np.sum(((linFunc(zx, b, -m) - zm)) ** 2)\n",
- " pparm = dict(throw_nan=False, pedantic=False,\n",
- " print_level=0)\n",
+ " \n",
+ " pparm = dict(throw_nan=False, pedantic=False, print_level=0)\n",
" pparm[\"m\"] = 0\n",
- " pparm[\"limit_m\"] = (0.9, 1)\n",
+ " if gain == 1:\n",
+ " pparm[\"limit_m\"] = (0.9, 1)\n",
+ " elif gain == 64:\n",
+ " pparm[\"limit_m\"] = (0, 1)\n",
" pparm[\"b\"] = np.nanmean(zm)\n",
"\n",
" mini = Minuit(linWrap, **pparm)\n",
" fmin = mini.migrad()\n",
- "\n",
- "\n",
- " if fmin[0]['is_valid'] and fmin[0]['has_covariance'] and (\n",
- " min_data is None or zx.size > min_data):\n",
- "\n",
- " ferr = mini.minos()\n",
+ " \n",
+ " if fmin[0]['is_valid'] and fmin[0]['has_covariance'] and (min_data is None or zx.size > min_data):\n",
+ " \n",
+ " ferr = mini.minos(maxcall = 10000)\n",
" res = mini.values\n",
"\n",
" ms[i] = res[\"m\"]\n",
@@ -754,36 +823,37 @@
" mErr[i] = max(-ferr['m']['lower'], ferr['m']['upper'])\n",
" bErr[i] = max(-ferr['b']['lower'], ferr['b']['upper'])\n",
"\n",
- " \n",
" zx, uidx = np.unique(zx, return_index=True)\n",
" zm = zm[uidx]\n",
" \n",
" yt = linFunc(zx, res[\"b\"], res[\"m\"])\n",
- " \n",
- " sidx = np.in1d(rtest, zx)\n",
- " \n",
- " entry['residuals'] = residuals\n",
+ " \n",
+ " else:\n",
+ " counter += 1 \n",
" \n",
- " ctiErr = np.sqrt((1. / bs * mErr) ** 2 + (ms / bs ** 2 * bErr) ** 2)\n",
" cti = ms / bs\n",
+ " ctiErr = np.sqrt((1. / bs * mErr) ** 2 + (ms / bs ** 2 * bErr) ** 2)\n",
" \n",
- " relGain = bs/np.nanmedian(bs)\n",
+ " relGain = bs/np.nanmean(bs) # I am not sure why Steffen was using median instead of mean here!!\n",
+ " relGainErr = np.sqrt((1. / np.nanmean(bs) * bErr) ** 2 + (bs / np.nanmean(bs) ** 2 * np.std(bs)) ** 2)\n",
+ "\n",
" entry['cti'] = cti\n",
" entry['cti_err'] = ctiErr\n",
" entry['rel_gain'] = relGain\n",
- " \n",
- " "
+ " entry['rel_gain_err'] = relGainErr\n",
+ "\n",
+ "print(\"Fitting is performed for the entire detector, and CTI and relative gain are calculated.\")\n",
+ "print(f\"Total number of fits which did not converge: {counter}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": false
- },
+ "metadata": {},
"outputs": [],
"source": [
- "\n",
+ "display(Markdown('### CTI per Quadrant'))\n",
+ " \n",
"fig = plt.figure(figsize=(15., 15.))\n",
"grid = AxesGrid(fig, 111, # similar to subplot(111)\n",
" nrows_ncols=(2, 2), # creates 2x2 grid of axes\n",
@@ -799,10 +869,13 @@
" cti = entry[\"cti\"]\n",
" x = np.arange(cti.size)+co[2]\n",
" grid[i].scatter(cti, x, 1., color='k')\n",
- " grid[i].set_xlim(1e-4, 1e-3)\n",
+ " if gain == 1:\n",
+ " grid[i].set_xlim(1e-4, 1e-3)\n",
+ " elif gain == 64:\n",
+ " grid[i].set_xlim(1e-18, 1e-1)\n",
" grid[i].semilogx()\n",
" grid[i].set_xlabel(\"CTI\")\n",
- " grid[i].set_ylabel(\"Row coordinate\")\n",
+ " grid[i].set_ylabel(\"Row Coordinate\")\n",
" grid[i].annotate(key, (0.1, 1.02), xycoords=\"axes fraction\", color=\"red\", annotation_clip=False)\n",
" grid[i].set_ylim(co[3], co[2])\n",
" ai += 1"
@@ -811,11 +884,45 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": false
- },
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "display(Markdown('### Uncertainty in CTI per Quadrant'))\n",
+ " \n",
+ "fig = plt.figure(figsize=(15., 15.))\n",
+ "grid = AxesGrid(fig, 111, # similar to subplot(111)\n",
+ " nrows_ncols=(2, 2), # creates 2x2 grid of axes\n",
+ " axes_pad=0.3, # pad between axes in inch.\n",
+ " aspect=False,\n",
+ " )\n",
+ "\n",
+ "ais = [0, 2, 1, 3]\n",
+ "ai = 0\n",
+ "for key, entry in dstats.items():\n",
+ " i = ais[ai]\n",
+ " co = entry[\"coords\"]\n",
+ " cti_err = entry[\"cti_err\"]\n",
+ " x = np.arange(cti_err.size)+co[2]\n",
+ " grid[i].scatter(cti_err, x, 1., color='k')\n",
+ " if gain == 1:\n",
+ " grid[i].set_xlim(1e-11, 1e-8)\n",
+ " elif gain == 64:\n",
+ " grid[i].set_xlim(1e-10, 1e-5)\n",
+ " grid[i].semilogx()\n",
+ " grid[i].set_xlabel(\"Uncertainty in CTI\")\n",
+ " grid[i].set_ylabel(\"Row Coordinate\")\n",
+ " grid[i].annotate(key, (0.1, 1.02), xycoords=\"axes fraction\", color=\"red\", annotation_clip=False)\n",
+ " grid[i].set_ylim(co[3], co[2])\n",
+ " ai += 1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
"outputs": [],
"source": [
+ "display(Markdown('### Relative Gain per Quadrant'))\n",
"\n",
"fig = plt.figure(figsize=(15., 15.))\n",
"grid = AxesGrid(fig, 111, # similar to subplot(111)\n",
@@ -833,51 +940,201 @@
" x = np.arange(rgain.size)+co[2]\n",
" grid[i].scatter(rgain, x, 1., color='k')\n",
" grid[i].set_xlim(0.5, 1.5)\n",
- " grid[i].set_xlabel(\"Relative gain\")\n",
- " grid[i].set_ylabel(\"Row coordinate\")\n",
+ " grid[i].set_xlabel(\"Relative Gain\")\n",
+ " grid[i].set_ylabel(\"Row Coordinate\")\n",
+ " grid[i].annotate(key, (0.1, 1.02), xycoords=\"axes fraction\", color=\"red\", annotation_clip=False)\n",
+ " grid[i].set_ylim(co[3], co[2])\n",
+ " ai += 1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "display(Markdown('### Uncertainty in Relative Gain per Quadrant'))\n",
+ "\n",
+ "fig = plt.figure(figsize=(15., 15.))\n",
+ "grid = AxesGrid(fig, 111, # similar to subplot(111)\n",
+ " nrows_ncols=(2, 2), # creates 2x2 grid of axes\n",
+ " axes_pad=0.3, # pad between axes in inch.\n",
+ " aspect=False,\n",
+ " )\n",
+ "\n",
+ "ais = [0, 2, 1, 3]\n",
+ "ai = 0\n",
+ "for key, entry in dstats.items():\n",
+ " i = ais[ai]\n",
+ " co = entry[\"coords\"]\n",
+ " rgain_err = entry[\"rel_gain_err\"]\n",
+ " x = np.arange(rgain_err.size)+co[2]\n",
+ " grid[i].scatter(rgain_err, x, 1., color='k')\n",
+ " grid[i].set_xlim(0, 0.2)\n",
+ " grid[i].set_xlabel(\"Relative Gain Uncertainty\")\n",
+ " grid[i].set_ylabel(\"Row Coordinate\")\n",
" grid[i].annotate(key, (0.1, 1.02), xycoords=\"axes fraction\", color=\"red\", annotation_clip=False)\n",
" grid[i].set_ylim(co[3], co[2])\n",
" ai += 1"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Sum of the total number of bad columns in each quadrant should be equal to the value of \"counter\"\n",
+ "\n",
+ "didnt_converge = [] # list of columns in each quadrant for which the fit did not converge. Here, we assume the \n",
+ " # number of columns in each quadrant always starts from 0 and ends at 512\n",
+ " \n",
+ "print(\"Rows for which fits did not converge:\\n\")\n",
+ "\n",
+ "for key, entry in dstats.items():\n",
+ " rel_gain = entry[\"rel_gain\"] \n",
+ " print(\"{}:\".format(key))\n",
+ " for index, element in enumerate(rel_gain):\n",
+ " if element == 0:\n",
+ " didnt_converge.append(index)\n",
+ " print(\"Total number: {}\".format(len(didnt_converge)))\n",
+ " print(\"These rows are as follows:\")\n",
+ " print(didnt_converge)\n",
+ " print(\"\\n\")\n",
+ " didnt_converge = []"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
"metadata": {
- "collapsed": false
+ "scrolled": false
},
"outputs": [],
"source": [
+ "display(Markdown('### CTE, CTI, and Relative Gain Maps'))\n",
+ "\n",
+ "ctemap = np.zeros((sensorSize[0], sensorSize[1]))\n",
+ "ctimap = np.zeros((sensorSize[0], sensorSize[1]))\n",
"gainmap = np.zeros((sensorSize[0], sensorSize[1]))\n",
- "for key, entry in dstats.items():\n",
+ "productmap = np.zeros((sensorSize[0], sensorSize[1]))\n",
"\n",
+ "for key, entry in dstats.items():\n",
" co = entry[\"coords\"]\n",
- " rel_gain = entry[\"rel_gain\"]\n",
+ " rel_gain = entry[\"rel_gain\"] \n",
" cti = entry[\"cti\"]\n",
" if entry[\"col_dir\"] == 1:\n",
- " quadgain = rel_gain[:, None]*(1-cti[:, None])**np.repeat(np.arange(512)[None, :], 512, axis=0)\n",
+ " # plus sign is so that the charge transfer inefficiency is 100% at the center of the CCD, which is the \n",
+ " # farthest distance away from the readout. If using a minus sign, you can set CTE = 100% at the readout \n",
+ " # => most efficient on the readout and less efficient towards the center of the chip:\n",
+ " quadcte = (1-cti[:, None])**np.repeat(np.arange(512)[None, :], 512, axis=0)\n",
" else:\n",
- " quadgain = rel_gain[:, None]*(1-cti[:, None])**np.repeat(np.arange(512, 0, -1)[None, :], 512, axis=0)\n",
- " gainmap[co[2]:co[3], co[0]:co[1]] = quadgain\n",
- "\n",
- "fig = xana.heatmapPlot(gainmap, figsize=(8,8), lut_label=\"Relative gain\")"
+ " quadcte = (1-cti[:, None])**np.repeat(np.arange(512, 0, -1)[None, :], 512, axis=0)\n",
+ " \n",
+ " quadcti = cti[:, None] \n",
+ " quadgain = rel_gain[:, None]\n",
+ " quadproduct = quadgain*quadcte\n",
+ " \n",
+ " ctemap[co[2]:co[3], co[0]:co[1]] = quadcte\n",
+ " ctimap[co[2]:co[3], co[0]:co[1]] = quadcti\n",
+ " gainmap[co[2]:co[3], co[0]:co[1]] = quadgain \n",
+ " productmap[co[2]:co[3], co[0]:co[1]] = quadproduct\n",
+ " \n",
+ "fig = xana.heatmapPlot(ctemap, figsize=(8, 8), x_label='Columns', y_label='Rows', \n",
+ " lut_label='Charge Transfer Efficiency', \n",
+ " aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), cmap = 'viridis',\n",
+ " panel_x_label='Along Rows', panel_y_label='Along Columns', \n",
+ " title = 'CTE Map for pnCCD (Gain = {})'.format(gain))\n",
+ " \n",
+ "fig = xana.heatmapPlot(ctimap, figsize=(8, 8), x_label='Columns', y_label='Rows', \n",
+ " lut_label='Charge Transfer Inefficiency',\n",
+ " norm=matplotlib.colors.LogNorm(vmin=np.nanmin(ctimap), vmax=np.nanmax(ctimap)),\n",
+ " aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), \n",
+ " panel_x_label='Along Rows', panel_y_label='Along Columns', \n",
+ " title = 'CTI Map for pnCCD (Gain = {})'.format(gain))\n",
+ "\n",
+ "fig = xana.heatmapPlot(gainmap, figsize=(8, 8), x_label='Columns', y_label='Rows', \n",
+ " lut_label='Relative Gain', vmin=0.8, vmax=1.2,\n",
+ " aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x),\n",
+ " panel_top_low_lim = 0.9, panel_top_high_lim = 1.1,\n",
+ " panel_x_label='Along Rows', panel_y_label='Along Columns', \n",
+ " title = 'Relativ Gain Map for pnCCD (Gain = {})'.format(gain))\n",
+ "\n",
+ "fig = xana.heatmapPlot(productmap, figsize=(8, 8), x_label='Columns', y_label='Rows', \n",
+ " lut_label='Relative Gain*CTI', \n",
+ " aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=0.8, vmax=1.2, \n",
+ " panel_x_label='Along Rows', panel_y_label='Along Columns', \n",
+ " panel_top_low_lim = 0.5, panel_top_high_lim = 1.5, panel_side_low_lim = 0.5, \n",
+ " panel_side_high_lim = 1.5, title = 'Relative Gain*CTE Map for pnCCD (Gain = {})'\n",
+ " .format(gain))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Sending the CTI and Relative Gain Maps to the Calibration Database"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": true
- },
+ "metadata": {},
"outputs": [],
- "source": []
+ "source": [
+ "constant_maps = {\n",
+ " 'RelativeGain' : gainmap,\n",
+ " 'CTI' : ctimap}\n",
+ "\n",
+ "for cname in constant_maps.keys():\n",
+ " metadata = ConstantMetaData()\n",
+ " det = Constants.CCD(DetectorTypes.pnCCD)\n",
+ " const = getattr(det, cname)()\n",
+ " const.data = constant_maps[cname].data\n",
+ " \n",
+ " metadata.calibration_constant = const\n",
+ " \n",
+ " # setting the operating condition\n",
+ " condition = Conditions.Illuminated.CCD(bias_voltage=bias_voltage,\n",
+ " photon_energy=photon_energy,\n",
+ " integration_time=integration_time,\n",
+ " gain_setting=gain,\n",
+ " temperature=fix_temperature,\n",
+ " pixels_x=pixels_x,\n",
+ " pixels_y=pixels_y)\n",
+ "\n",
+ " device = Detectors.PnCCD1\n",
+ " metadata.detector_condition = condition\n",
+ " \n",
+ " # Specifying 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",
+ " \n",
+ " metadata.calibration_constant_version.raw_data_location = file_loc\n",
+ " \n",
+ " if db_output:\n",
+ " try:\n",
+ " metadata.send(cal_db_interface, timeout=cal_db_timeout)\n",
+ " print(f\"Inject {cname} constants from {metadata.calibration_constant_version.begin_at}\")\n",
+ " except Exception as e:\n",
+ " print(f\"Error: {e}\")\n",
+ " \n",
+ " if local_output:\n",
+ " save_const_to_h5(metadata, out_folder)\n",
+ " print(f\"Calibration constant {cname} is stored to {out_folder}.\")\n",
+ " \n",
+ "print(\"\\nGenerated constants with conditions:\\n\")\n",
+ "print(f\"• bias_voltage: {bias_voltage}\\n• photon_energy: {photon_energy}\\n\"\n",
+ " f\"• integration_time: {integration_time}\\n• gain_setting: {gain}\\n\"\n",
+ " f\"• temperature: {fix_temperature}\\n• creation_time: {creation_time}\\n\")"
+ ]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "collapsed": true
- },
+ "metadata": {},
"outputs": [],
"source": []
}
@@ -898,7 +1155,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.6.6"
+ "version": "3.6.7"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
diff --git a/notebooks/pnCCD/Correct_pnCCD_NBC.ipynb b/notebooks/pnCCD/Correct_pnCCD_NBC.ipynb
index b4ddf1567..8723b6ef5 100644
--- a/notebooks/pnCCD/Correct_pnCCD_NBC.ipynb
+++ b/notebooks/pnCCD/Correct_pnCCD_NBC.ipynb
@@ -4,11 +4,11 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# pnCCD Data Correction ##\n",
+ "# pnCCD Data Correction #\n",
"\n",
- "Authors: DET Group, Version 1.0\n",
+ "Authors: DET Group, Modified by Kiana Setoodehnia on March 2020 - Version 2.0\n",
"\n",
- "The following notebook provides correction of images acquired with the pnCCD."
+ "The following notebook provides offset, relative gain, common mode, split events, and pattern classification corrections of images acquired with the pnCCD. This notebook *does not* yet correct for charge transfer inefficiency."
]
},
{
@@ -22,44 +22,69 @@
},
"outputs": [],
"source": [
- "cluster_profile = \"noDB\" # ipcluster profile to use\n",
- "in_folder = \"/gpfs/exfel/exp/SQS/201921/p002430/raw/\" # input folder, required\n",
- "out_folder = '/gpfs/exfel/data/scratch/karnem/test/' # output folder, required\n",
+ "in_folder = \"/gpfs/exfel/exp/SQS/201930/p900075/raw\" # input folder\n",
+ "out_folder = '/gpfs/exfel/data/scratch/ahmedk/test/pnccd' # output folder\n",
+ "run = 365 # which run to read data from\n",
"sequences = [-1] # sequences to correct, set to -1 for all, range allowed\n",
- "run = 746 # which run to read data from, required\n",
"\n",
"karabo_da = 'PNCCD01' # data aggregators\n",
"karabo_id = \"SQS_NQS_PNCCD1MP\" # karabo prefix of PNCCD devices\n",
"receiver_id = \"PNCCD_FMT-0\" # inset for receiver devices\n",
- "path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data\n",
- "h5path = '/INSTRUMENT/{}/CAL/{}:output/data/' # path in the HDF5 file the data is at\n",
- "\n",
- "use_dir_creation_date = True\n",
- "cal_db_interface = \"tcp://max-exfl016:8015\" # calibration DB interface to use\n",
- "cal_db_timeout = 300000000 # timeout on caldb requests\n",
+ "path_template = 'RAW-R{:04d}-PNCCD01-S{{:05d}}.h5' # the template to use to access data\n",
+ "path_template_seqs = \"{}/r{:04d}/*PNCCD01-S*.h5\"\n",
+ "h5path = '/INSTRUMENT/{}/CAL/{}:output/data/' # path to data in the HDF5 file \n",
"\n",
+ "overwrite = True # keep this as True to not overwrite the output \n",
+ "use_dir_creation_date = True # required to obtain creation time of the run\n",
"number_dark_frames = 0 # number of images to be used, if set to 0 all available images are used\n",
- "sigma_noise = 10. # Pixel exceeding 'sigmaNoise' * noise value in that pixel will be masked\n",
- "multi_iteration = False # use multiple iterations\n",
- "fix_temperature = 233.\n",
- "gain = 0\n",
- "bias_voltage = 300\n",
- "integration_time = 70\n",
+ "chunk_size_idim = 1 # H5 chunking size of output data\n",
+ "cluster_profile = \"noDB\" # ipcluster profile to use\n",
+ "\n",
+ "cpuCores = 8\n",
+ "commonModeBlockSize = [512, 512] # size of the detector in pixels for common mode calculations\n",
+ "commonModeAxis = 0 # axis along which common mode will be calculated, 0 = row, and 1 = column \n",
"split_evt_primary_threshold = 5. # primary threshold for split event classification in terms of n sigma noise\n",
"split_evt_secondary_threshold = 3. # secondary threshold for split event classification in terms of n sigma noise\n",
"split_evt_mip_threshold = 1000. # MIP threshold for event classification\n",
- "chunk_size_idim = 1 # H5 chunking size of output data\n",
- "overwrite = True\n",
- "do_pattern_classification = True # classify split events\n",
"sequences_per_node = 1\n",
"limit_images = 0\n",
- "time_offset_days = 0\n",
- "photon_energy_gain_map = 2. # energy in keV\n",
- "cpuCores = 8\n",
+ "seq_num = 0 # sequence number for which the last plot at the end of the notebook is plotted\n",
+ "\n",
+ "# pnCCD parameters:\n",
+ "fix_temperature = 233.\n",
+ "gain = 1\n",
+ "bias_voltage = 300\n",
+ "integration_time = 70\n",
+ "photon_energy = 1.6 # Al fluorescence in keV\n",
+ "\n",
+ "cal_db_interface = \"tcp://max-exfl016:8015\" # calibration DB interface to use\n",
+ "cal_db_timeout = 300000 # timeout on caldb requests\n",
+ "creation_time = \"\" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00\n",
+ "\n",
+ "only_offset = False # Only, apply offset.\n",
+ "common_mode = True # Apply common mode correction\n",
+ "relgain = True # Apply relative gain correction\n",
+ "cti = False # Apply charge transfer inefficiency correction (not implemented, yet)\n",
+ "do_pattern_classification = True # classify split events"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Here the herarichy and dependability for correction booleans are defined \n",
+ "corr_bools = {}\n",
"\n",
- "def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da):\n",
- " from xfel_calibrate.calibrate import balance_sequences as bs\n",
- " return bs(in_folder, run, sequences, sequences_per_node, karabo_da)"
+ "corr_bools[\"only_offset\"] = only_offset\n",
+ "\n",
+ "# Dont apply any corrections if only_offset is requested.\n",
+ "if not only_offset:\n",
+ " corr_bools[\"relgain\"] = relgain\n",
+ " corr_bools[\"cti\"] = cti\n",
+ " corr_bools[\"common_mode\"] = common_mode\n",
+ " corr_bools[\"pattern_class\"] = do_pattern_classification"
]
},
{
@@ -73,136 +98,144 @@
},
"outputs": [],
"source": [
- "import XFELDetAna.xfelprofiler as xprof\n",
+ "import warnings\n",
+ "warnings.filterwarnings('ignore')\n",
+ "import h5py\n",
+ "import time\n",
+ "import copy\n",
+ "import os\n",
+ "import glob\n",
+ "import gc\n",
+ "import datetime\n",
+ "from datetime import timedelta\n",
+ "from prettytable import PrettyTable\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "from iminuit import Minuit\n",
+ "from IPython.display import display, Markdown\n",
"\n",
+ "import XFELDetAna.xfelprofiler as xprof\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 glob\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",
+ "from cal_tools.tools import get_dir_creation_date, get_random_db_interface\n",
"\n",
"%matplotlib inline\n",
"\n",
"if sequences[0] == -1:\n",
- " sequences = None\n",
- "\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)))])"
+ " sequences = None"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2018-12-06T15:54:23.679069Z",
- "start_time": "2018-12-06T15:54:23.662821Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "\n",
- "x = 1024 # rows of the FastCCD to analyze in FS mode \n",
- "y = 1024 # columns of the FastCCD to analyze in FS mode \n",
- "\n",
+ "# Each xcal.HistogramCalculator requires a total number of bins and a binning range. We define these using a \n",
+ "# dictionary:\n",
+ "\n",
+ "# For all xcal histograms:\n",
+ "if gain == 1:\n",
+ " Hist_Bin_Dict = {\n",
+ " \"bins\": 70000, # number of bins \n",
+ " \"bin_range\": [0, 70000]\n",
+ " }\n",
+ "\n",
+ " # For the numpy histograms on the last cell of the notebook:\n",
+ " Event_Bin_Dict = {\n",
+ " \"event_bins\": 1000, # number of bins \n",
+ " \"b_range\": [0, 50000] # bin range \n",
+ " }\n",
" \n",
- "ped_dir = \"{}/r{:04d}\".format(in_folder, run)\n",
- "fp_name = path_template.format(run, karabo_da)\n",
- "h5path = h5path.format(karabo_id, receiver_id)\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()))\n"
+ "elif gain == 64:\n",
+ " # For all xcal histograms:\n",
+ " Hist_Bin_Dict = {\n",
+ " \"bins\": 25000, # number of bins \n",
+ " \"bin_range\": [0, 25000] \n",
+ " }\n",
+ " # For the numpy histograms on the last cell of the notebook:\n",
+ " Event_Bin_Dict = {\n",
+ " \"event_bins\": 1000, # number of bins \n",
+ " \"b_range\": [0, 3000] # bin range \n",
+ " }\n",
+ " \n",
+ "bins = Hist_Bin_Dict[\"bins\"]\n",
+ "bin_range = Hist_Bin_Dict[\"bin_range\"]\n",
+ "event_bins = Event_Bin_Dict[\"event_bins\"]\n",
+ "b_range = Event_Bin_Dict[\"b_range\"]\n",
+ "\n",
+ "# On the singles spectrum (uploaded in the middle of this notebook), the ADU values correspoding to the boundaries\n",
+ "# of the first peak region are used as cti_limit_low and cti_limit_high:\n",
+ "\n",
+ "if gain == 1:\n",
+ " cti_limit_low = 3000 # lower limit of cti\n",
+ " cti_limit_high = 10000 # higher limit of cti\n",
+ "elif gain == 64:\n",
+ " cti_limit_low = 50\n",
+ " cti_limit_high = 170"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2018-12-06T15:54:23.913269Z",
- "start_time": "2018-12-06T15:54:23.868910Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
+ "# Calibration Database Settings, and Some Initial Run Parameters & Paths:\n",
"\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 = [512, 512]\n",
- "commonModeAxisR = 'row'#Axis along which common mode will be calculated\n",
- "run_parallel = True\n",
- "profile = False\n",
+ "display(Markdown('### Initial Settings and Paths'))\n",
+ "pixels_x = 1024 # rows of pnCCD in pixels \n",
+ "pixels_y = 1024 # columns of pnCCD in pixels\n",
+ "print(f\"pnCCD size is: {pixels_x}x{pixels_y} pixels.\")\n",
+ "print(f'Calibration database interface selected: {cal_db_interface}')\n",
+ "\n",
+ "proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]\n",
+ "print(f'Proposal: {proposal}, Run: {run}')\n",
+ "\n",
+ "# Paths to the data:\n",
+ "ped_dir = \"{}/r{:04d}\".format(in_folder, run)\n",
+ "fp_name = path_template.format(run, karabo_da)\n",
+ "fp_path = '{}/{}'.format(ped_dir, fp_name)\n",
+ "h5path = h5path.format(karabo_id, receiver_id)\n",
+ "print(\"HDF5 path to data: {}\\n\".format(h5path))\n",
+ "\n",
+ "# Run's creation time:\n",
+ "if creation_time:\n",
+ " try:\n",
+ " creation_time = datetime.datetime.strptime(creation_time, '%Y-%m-%d %H:%M:%S.%f')\n",
+ " except Exception as e:\n",
+ " print(f\"creation_time value error: {e}.\" \n",
+ " \"Use same format as YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00/n\")\n",
+ " creation_time = None\n",
+ " print(\"Given creation time wont be used.\")\n",
+ "else:\n",
+ " creation_time = None\n",
+ "\n",
+ "if not creation_time and use_dir_creation_date:\n",
+ " creation_time = get_dir_creation_date(in_folder, run)\n",
+ "\n",
+ "\n",
+ "print(f\"Creation time: {creation_time}\")\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"
- }
- },
- "outputs": [],
- "source": [
+ "# Reading all sequences of the run:\n",
"file_list = []\n",
"total_sequences = 0\n",
"fsequences = []\n",
"\n",
"if sequences is None:\n",
" file_list = glob.glob(fp_path.format(0).replace('00000', '*'))\n",
+ " file_list = sorted(file_list, key = lambda x: (len (x), x))\n",
" total_sequences = len(file_list)\n",
" fsequences = range(total_sequences)\n",
"else:\n",
@@ -213,7 +246,20 @@
" total_sequences += 1\n",
" fsequences.append(seq)\n",
"\n",
- "sequences = fsequences"
+ "sequences = fsequences \n",
+ "print(f\"This run has a total number of {total_sequences} sequences.\\n\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "display(Markdown('### List of Files to be Processed'))\n",
+ "print(\"Reading data from the following files:\")\n",
+ "for i in range(len(file_list)):\n",
+ " print(file_list[i])"
]
},
{
@@ -221,30 +267,18 @@
"execution_count": null,
"metadata": {
"ExecuteTime": {
- "end_time": "2018-12-06T18:43:39.776018Z",
- "start_time": "2018-12-06T18:43:39.759185Z"
+ "end_time": "2018-12-06T15:54:23.913269Z",
+ "start_time": "2018-12-06T15:54:23.868910Z"
}
},
"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",
+ "# Sensor size and block size definitions (important for common mode and other corrections):\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."
+ "sensorSize = [pixels_x, pixels_y]\n",
+ "blockSize = [sensorSize[0]//2, sensorSize[1]//2] # sensor area will be analysed according to blockSize\n",
+ "xcal.defaultBlockSize = blockSize # for xcal.HistogramCalculators \n",
+ "memoryCells = 1 # pnCCD has 1 memory cell"
]
},
{
@@ -252,247 +286,203 @@
"execution_count": null,
"metadata": {
"ExecuteTime": {
- "end_time": "2018-12-06T15:54:28.254544Z",
- "start_time": "2018-12-06T15:54:24.709521Z"
+ "end_time": "2018-12-06T15:54:23.913269Z",
+ "start_time": "2018-12-06T15:54:23.868910Z"
}
},
"outputs": [],
"source": [
- "offsetMap = None\n",
- "offset_temperature = None\n",
- "badPixelMap = None\n",
- "noiseMap = None\n",
- "\n",
- "## offset\n",
- "metadata = ConstantMetaData()\n",
- "offset = Constants.CCD(DetectorTypes.pnCCD).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=gain,\n",
- " temperature=fix_temperature,\n",
- " pixels_x=1024,\n",
- " pixels_y=1024)\n",
- "\n",
- "\n",
- "device = Detectors.PnCCD1\n",
- "\n",
- "\n",
- "metadata.detector_condition = condition\n",
- "\n",
- "\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",
- " 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 = 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 dark images for offset calculation \" +\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.pnCCD).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=gain,\n",
- " temperature=fix_temperature,\n",
- " pixels_x=1024,\n",
- " pixels_y=1024)\n",
- "\n",
- "\n",
- "device = Detectors.PnCCD1\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",
- " 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",
- "noiseMap = noise.data\n",
- "\n",
- "## bad pixels \n",
- "\n",
- "metadata = ConstantMetaData()\n",
- "bpix = Constants.CCD(DetectorTypes.pnCCD).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=gain,\n",
- " temperature=fix_temperature,\n",
- " pixels_x=1024,\n",
- " pixels_y=1024)\n",
- "device = Detectors.PnCCD1\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, timeout=cal_db_timeout)\n",
- "\n",
- " \n",
- "badPixelMap = bpix.data"
+ "# Output Folder Creation:\n",
+ "os.makedirs(out_folder, exist_ok=True)\n",
+ "if not overwrite:\n",
+ " raise AttributeError(\"Output path exists! Exiting\") "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Loading cti and relative gain values"
+ "As a first step, dark constants have to be retrieved from the calibration database"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2018-12-06T15:54:28.343869Z",
- "start_time": "2018-12-06T15:54:28.271344Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "## relative gain\n",
- "\n",
- "metadata = ConstantMetaData()\n",
- "relgain = Constants.CCD(DetectorTypes.pnCCD).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=gain,\n",
- " temperature=fix_temperature,\n",
- " pixels_x=1024,\n",
- " pixels_y=1024, photon_energy=photon_energy_gain_map)\n",
- "device = Detectors.PnCCD1\n",
- "\n",
+ "def get_dark(db_parms, bias_voltage, integration_time, gain, fix_temperature, creation_time, run):\n",
+ "# This function is to retrieve the dark constants associated with the run of interest:\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 = np.ones_like(offsetMap) #relgain.data\n"
+ " cal_db_interface, cal_db_timeout = db_parms\n",
+ " cal_db_interface = get_random_db_interface(cal_db_interface)\n",
+ " print(f'Calibration database interface: {cal_db_interface}')\n",
+ " \n",
+ " try: \n",
+ " Offset = None\n",
+ " Noise = None\n",
+ " BadPixelsDark = None\n",
+ "\n",
+ " constants = {} \n",
+ " constants['Offset'] = Offset\n",
+ " constants['Noise'] = Noise\n",
+ " constants['BadPixelsDark'] = BadPixelsDark\n",
+ "\n",
+ " for const in constants.keys():\n",
+ " metadata = ConstantMetaData()\n",
+ " dconst = getattr(Constants.CCD(DetectorTypes.pnCCD), const)()\n",
+ " dconst.data = constants[const]\n",
+ " metadata.calibration_constant = dconst\n",
+ "\n",
+ " condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,\n",
+ " integration_time=integration_time,\n",
+ " gain_setting=gain,\n",
+ " temperature=fix_temperature,\n",
+ " pixels_x=pixels_x,\n",
+ " pixels_y=pixels_y)\n",
+ "\n",
+ " device = Detectors.PnCCD1\n",
+ " metadata.detector_condition = condition\n",
+ "\n",
+ " metadata.calibration_constant_version = Versions.Timespan(device=device, start=creation_time)\n",
+ " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=cal_db_timeout)\n",
+ "\n",
+ " constants[const] = dconst.data \n",
+ " except Exception as e:\n",
+ " print(f\"Error: {e}\")\n",
+ " for parm in condition.parameters:\n",
+ " if parm.name == \"Sensor Temperature\":\n",
+ " print(f\"For run {run}: dark images for offset calculation \" \n",
+ " f\"were taken at temperature: {parm.value:0.2f} \"\n",
+ " f\"K @ {metadata.calibration_constant_version.begin_at}\") \n",
+ " return constants"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
- "ExecuteTime": {
- "end_time": "2018-12-06T15:54:28.771629Z",
- "start_time": "2018-12-06T15:54:28.346051Z"
- }
+ "scrolled": false
},
"outputs": [],
"source": [
- "#************************Calculators************************#\n",
- "\n",
- "cmCorrection = xcal.CommonModeCorrection([x, y],\n",
- " commonModeBlockSize,\n",
- " commonModeAxisR,\n",
- " parallel=False, dType=np.float32, stride=1,\n",
- " noiseMap=noiseMap.astype(np.float32), minFrac=0)\n",
- "# cmCorrection.debug()\n",
- "\n",
- "patternClassifierLH = xcal.PatternClassifier([x, y//2],\n",
- " noiseMap[:, :y//2],\n",
- " split_evt_primary_threshold,\n",
- " split_evt_secondary_threshold,\n",
- " split_evt_mip_threshold,\n",
- " tagFirstSingles=3,\n",
- " nCells=memoryCells,\n",
- " cores=cpuCores,\n",
- " allowElongated=False,\n",
- " blockSize=[x, y//2],\n",
- " runParallel=True)\n",
- "\n",
- "\n",
- "patternClassifierUH = xcal.PatternClassifier([x, y//2],\n",
- " noiseMap[:, y//2:],\n",
- " split_evt_primary_threshold,\n",
- " split_evt_secondary_threshold,\n",
- " split_evt_mip_threshold,\n",
- " tagFirstSingles=4,\n",
- " nCells=memoryCells,\n",
- " cores=cpuCores,\n",
- " allowElongated=False,\n",
- " blockSize=[x, y//2],\n",
- " runParallel=True)"
+ "display(Markdown('### Dark Data Retrieval'))\n",
+ "\n",
+ "db_parms = cal_db_interface, cal_db_timeout\n",
+ "\n",
+ "constants = get_dark(db_parms, bias_voltage, integration_time, gain, fix_temperature, creation_time, run)\n",
+ "\n",
+ "fig = xana.heatmapPlot(constants[\"Offset\"][:,:,0], x_label='Columns', y_label='Rows', lut_label='Pedestal (ADU)', aspect=1, \n",
+ " x_range=(0, pixels_y), y_range=(0, pixels_x), vmax=16000, \n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " title = 'Dark Offset Map')\n",
+ "\n",
+ "fig = xana.heatmapPlot(constants[\"Noise\"][:,:,0], x_label='Columns', y_label='Rows', \n",
+ " lut_label='Common Mode Corrected Noise (ADU)', \n",
+ " aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), # vmin=-50, vmax=70000, \n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " title = 'Dark Noise Map')\n",
+ "\n",
+ "fig = xana.heatmapPlot(np.log2(constants[\"BadPixelsDark\"][:,:,0]), x_label='Columns', y_label='Rows', \n",
+ " lut_label='Bad Pixel Value (ADU)', \n",
+ " aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), # vmin=-50, vmax=70000, \n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " title = 'Dark Bad Pixels Map')"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2018-12-06T16:08:51.886343Z",
- "start_time": "2018-12-06T16:08:51.842837Z"
- }
- },
+ "metadata": {},
"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)"
+ "if corr_bools.get('relgain'):\n",
+ " display(Markdown('We will now retrieve the relative gain map from the calibration database'))\n",
+ " metadata = ConstantMetaData()\n",
+ " relgain = Constants.CCD(DetectorTypes.pnCCD).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=gain,\n",
+ " temperature=fix_temperature,\n",
+ " pixels_x=pixels_x,\n",
+ " pixels_y=pixels_y, \n",
+ " photon_energy=photon_energy)\n",
+ "\n",
+ " device = Detectors.PnCCD1\n",
+ " metadata.detector_condition = condition\n",
+ "\n",
+ "\n",
+ " # specify the a version for this constant\n",
+ " metadata.calibration_constant_version = Versions.Now(device=device)\n",
+ " metadata.retrieve(cal_db_interface, timeout=cal_db_timeout) \n",
+ " constants[\"RelativeGain\"] = relgain.data\n",
+ " print(\"Retrieved RelativeGain constant with creation time:\"\n",
+ " f\"{metadata.calibration_constant_version.begin_at}\")\n",
+ " display(Markdown('### Relative Gain Map Retrieval'))\n",
+ " fig = xana.heatmapPlot(constants[\"RelativeGain\"], figsize=(8, 8), x_label='Columns', y_label='Rows', lut_label='Relative Gain', \n",
+ " aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=0.8, vmax=1.2, \n",
+ " panel_x_label='Along Rows', panel_y_label='Along Columns', \n",
+ " panel_top_low_lim = 0.5, panel_top_high_lim = 1.5, panel_side_low_lim = 0.5, \n",
+ " panel_side_high_lim = 1.5, \n",
+ " title = f'1st Relative Gain Map for pnCCD (Gain = {gain})')"
]
},
{
- "cell_type": "markdown",
+ "cell_type": "code",
+ "execution_count": null,
"metadata": {},
+ "outputs": [],
"source": [
- "Applying corrections"
+ "#************************ Calculators ************************#\n",
+ "\n",
+ "if corr_bools.get('common_mode'):\n",
+ " # Common Mode Correction Calculator:\n",
+ " cmCorrection = xcal.CommonModeCorrection([pixels_x, pixels_y],\n",
+ " commonModeBlockSize,\n",
+ " commonModeAxis,\n",
+ " parallel=False, dType=np.float32, stride=1,\n",
+ " noiseMap=constants[\"Noise\"].astype(np.float32), minFrac=0)\n",
+ " cmCorrection.debug()\n",
+ "\n",
+ "if corr_bools.get('pattern_class'):\n",
+ " # Pattern Classifier Calculator:\n",
+ " # Left Hemisphere:\n",
+ " patternClassifierLH = xcal.PatternClassifier([pixels_x, pixels_y//2],\n",
+ " constants[\"Noise\"][:, :pixels_y//2],\n",
+ " split_evt_primary_threshold,\n",
+ " split_evt_secondary_threshold,\n",
+ " split_evt_mip_threshold,\n",
+ " tagFirstSingles=3, # track along y-axis, left to right (see \n",
+ " nCells=memoryCells, # split_event.py file in pydetlib/lib/src/\n",
+ " cores=cpuCores, # XFELDetAna/algorithms)\n",
+ " allowElongated=False,\n",
+ " blockSize=[pixels_x, pixels_y//2],\n",
+ " runParallel=True)\n",
+ "\n",
+ " # Right Hemisphere:\n",
+ " patternClassifierRH = xcal.PatternClassifier([pixels_x, pixels_y//2],\n",
+ " constants[\"Noise\"][:, pixels_y//2:],\n",
+ " split_evt_primary_threshold,\n",
+ " split_evt_secondary_threshold,\n",
+ " split_evt_mip_threshold,\n",
+ " tagFirstSingles=4, # track along y-axis, right to left\n",
+ " nCells=memoryCells,\n",
+ " cores=cpuCores,\n",
+ " allowElongated=False,\n",
+ " blockSize=[pixels_x, pixels_y//2],\n",
+ " runParallel=True)\n",
+ "\n",
+ " patternClassifierLH._imagesPerChunk = 500\n",
+ " patternClassifierRH._imagesPerChunk = 500\n",
+ " patternClassifierLH.debug()\n",
+ " patternClassifierRH.debug()\n",
+ "\n",
+ " # Setting bad pixels:\n",
+ " patternClassifierLH.setBadPixelMask(constants[\"BadPixelsDark\"][:, :pixels_y//2] != 0)\n",
+ " patternClassifierRH.setBadPixelMask(constants[\"BadPixelsDark\"][:, pixels_y//2:] != 0)"
]
},
{
@@ -500,31 +490,70 @@
"execution_count": null,
"metadata": {
"ExecuteTime": {
- "end_time": "2018-12-06T16:08:52.441784Z",
- "start_time": "2018-12-06T16:08:52.437284Z"
+ "end_time": "2018-12-06T15:54:28.771629Z",
+ "start_time": "2018-12-06T15:54:28.346051Z"
}
},
"outputs": [],
"source": [
- "patternClassifierLH._imagesPerChunk = 500\n",
- "patternClassifierUH._imagesPerChunk = 500\n",
- "#patternClassifierLH.debug()\n",
- "#patternClassifierUH.debug()"
+ "#***************** Histogram Calculators ******************#\n",
+ "# Will contain uncorrected data:\n",
+ "histCalRaw = xcal.HistogramCalculator(sensorSize, \n",
+ " bins=bins, \n",
+ " range=bin_range,\n",
+ " nCells=memoryCells, \n",
+ " cores=cpuCores,\n",
+ " blockSize=blockSize) \n",
+ "histCalRaw.debug()\n",
+ "# Will contain offset corrected data:\n",
+ "histCalOffsetCor = xcal.HistogramCalculator(sensorSize, \n",
+ " bins=bins, \n",
+ " range=bin_range,\n",
+ " nCells=memoryCells, \n",
+ " cores=cpuCores,\n",
+ " blockSize=blockSize)\n",
+ "histCalOffsetCor.debug()\n",
+ "if corr_bools.get('common_mode'):\n",
+ " # Will contain common mode corrected data:\n",
+ " histCalCommonModeCor = xcal.HistogramCalculator(sensorSize, \n",
+ " bins=bins, \n",
+ " range=bin_range,\n",
+ " nCells=memoryCells, \n",
+ " cores=cpuCores,\n",
+ " blockSize=blockSize)\n",
+ " histCalCommonModeCor.debug()\n",
+ "# Will contain split events pattern data:\n",
+ "histCalPcorr = xcal.HistogramCalculator(sensorSize, \n",
+ " bins=bins, \n",
+ " range=bin_range,\n",
+ " nCells=memoryCells, \n",
+ " cores=cpuCores,\n",
+ " blockSize=blockSize)\n",
+ "histCalPcorr.debug()\n",
+ "# Will contain singles events data:\n",
+ "histCalPcorrS = xcal.HistogramCalculator(sensorSize, \n",
+ " bins=bins, \n",
+ " range=bin_range,\n",
+ " nCells=memoryCells, \n",
+ " cores=cpuCores,\n",
+ " blockSize=blockSize)\n",
+ "histCalPcorrS.debug()\n",
+ "if corr_bools.get('relgain'):\n",
+ " # Will contain gain corrected data:\n",
+ " histCalGainCor = xcal.HistogramCalculator(sensorSize, \n",
+ " bins=bins, \n",
+ " range=bin_range,\n",
+ " nCells=memoryCells, \n",
+ " cores=cpuCores,\n",
+ " blockSize=blockSize)\n",
+ " histCalGainCor.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"
- }
- },
- "outputs": [],
+ "cell_type": "markdown",
+ "metadata": {},
"source": [
- "histCalOffsetCor.debug()\n",
- "histCalPcorr.debug()\n"
+ "Applying offset and common mode corrections to the raw data"
]
},
{
@@ -539,12 +568,11 @@
"outputs": [],
"source": [
"def copy_and_sanitize_non_cal_data(infile, outfile, h5base):\n",
- " \n",
+ " '''This function reads the .h5 data and writes the corrected .h5 data.'''\n",
" if h5base.startswith(\"/\"):\n",
" h5base = h5base[1:]\n",
" dont_copy = ['image']\n",
- " dont_copy = [h5base+\"/{}\".format(do)\n",
- " for do in dont_copy]\n",
+ " dont_copy = [h5base+\"/{}\".format(do) for do in dont_copy]\n",
"\n",
" def visitor(k, item):\n",
" if k not in dont_copy:\n",
@@ -561,24 +589,25 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2018-12-06T16:10:55.917179Z",
- "start_time": "2018-12-06T16:09:01.603633Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "mean_im = None\n",
- "single_im = None\n",
+ "# Data corrections and event classifications happen here. Also, the corrected data are written to datasets:\n",
+ "uncor_mean_im = None\n",
+ "uncor_single_im = None\n",
+ "offset_mean_im = None\n",
+ "offset_single_im = None\n",
+ "cm_mean_im = None\n",
+ "cm_single_im = None\n",
+ "gain_mean_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",
+ "\n",
+ "offsetMap = np.squeeze(constants[\"Offset\"])\n",
+ "noiseMap = np.squeeze(constants[\"Noise\"])\n",
+ "badPixelMap = np.squeeze(constants[\"BadPixelsDark\"])\n",
+ "if corr_bools.get('relgain'):\n",
+ " relGain = constants[\"RelativeGain\"]\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",
@@ -586,17 +615,20 @@
"\n",
" data = None\n",
" noise = None\n",
+ " \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+\"/image\"][()]\n",
- " nzidx = np.count_nonzero(data, axis=(1,2))\n",
+ " nzidx = np.count_nonzero(data, axis=(1, 2))\n",
" data = data[nzidx != 0, ...]\n",
+ " \n",
" if limit_images > 0:\n",
" data = data[:limit_images,...]\n",
" oshape = data.shape\n",
" data = np.moveaxis(data, 0, 2)\n",
+ " \n",
" ddset = ofile.create_dataset(h5path+\"/pixels\",\n",
" oshape,\n",
" chunks=(chunk_size_idim, oshape[1], oshape[2]),\n",
@@ -606,178 +638,294 @@
" oshape,\n",
" chunks=(chunk_size_idim, oshape[1], oshape[2]),\n",
" dtype=np.uint32, compression=\"gzip\")\n",
+ "\n",
+ " if corr_bools.get('relgain'):\n",
+ " ddsetg = ofile.create_dataset(h5path+\"/gain\",\n",
+ " oshape,\n",
+ " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n",
+ " dtype=np.float32, compression=\"gzip\")\n",
" \n",
+ " data = data.astype(np.float32) \n",
" \n",
- " \n",
- " data = data.astype(np.float32) \n",
+ " # Creating maps for correction usage:\n",
" offset = np.repeat(offsetMap[...,None], data.shape[2], axis=2)\n",
- " rg = np.repeat(relGain[:,:,None], data.shape[2], axis=2)\n",
" noise = np.repeat(noiseMap[...,None], data.shape[2], axis=2)\n",
" bpix = np.repeat(badPixelMap[...,None], data.shape[2], axis=2)\n",
+ " if corr_bools.get('relgain'):\n",
+ " rg = np.repeat(relGain[:,:,None], data.shape[2], axis=2) # rg = relative gain\n",
" \n",
- " data -= offset\n",
- " data *= rg\n",
- " \n",
+ " # uncor_mean_im = averaged over all non-corrected images (in the first sequence only)\n",
+ " if uncor_mean_im is None:\n",
+ " uncor_mean_im = np.nanmean(data, axis=2) \n",
+ " uncor_single_im = data[...,0] # The non-corrected image corresponding to the first frame\n",
+ " \n",
+ " histCalRaw.fill(data) # filling histogram with raw uncorrected data\n",
" \n",
- " histCalOffsetCor.fill(data)\n",
+ " # masking data for bad pixels and equating the values to np.nan so that the pattern classifier \n",
+ " # ignores them:\n",
+ " data[bpix != 0] = np.nan\n",
+ " \n",
+ " data -= offset # offset correction \n",
+ " histCalOffsetCor.fill(data) # filling histogram with offset corrected data\n",
"\n",
- " \n",
" ddset[...] = np.moveaxis(data, 2, 0)\n",
" ddsetm[...] = np.moveaxis(bpix, 2, 0)\n",
" ofile.flush()\n",
- " \n",
- " if mean_im is None:\n",
- " mean_im = np.nanmean(data, axis=2)\n",
- " single_im = data[...,0]\n",
+ "\n",
+ " # mean_image = averaged over all offset corrected images (in the first sequence only)\n",
+ " if offset_mean_im is None:\n",
+ " offset_mean_im = np.nanmean(data, axis=2) \n",
+ " offset_single_im = data[...,0] # The offset corrected image corresponding to the first frame \n",
" \n",
- " if do_pattern_classification:\n",
+ " if corr_bools.get('relgain'):\n",
+ " data /= rg # relative gain correction \n",
+ " histCalGainCor.fill(data) # filling histogram with gain corrected data\n",
+ " ddsetg[...] = np.moveaxis(rg, 2, 0).astype(np.float32)\n",
+ "\n",
+ " # gain_mean_image = averaged over gain corrected images (in the first sequence only)\n",
+ " if gain_mean_im is None:\n",
+ " gain_mean_im = np.nanmean(data, axis=2) \n",
+ " gain_single_im = data[...,0] # The gain corrected image corresponding to the first frame\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",
+ " if corr_bools.get('pattern_class'):\n",
+ "\n",
+ " # cm: common mode, c: classifications, p: even patterns\n",
+ " if corr_bools.get('common_mode'):\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",
+ " 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",
+ " # The calculation of the cluster map:\n",
+ " patternClassifierLH._noisemap = noise[:, :pixels_x//2, :]\n",
+ " patternClassifierRH._noisemap = noise[:, pixels_x//2:, :]\n",
+ " if corr_bools.get('common_mode'):\n",
+ " data = cmCorrection.correct(data.astype(np.float32), # common mode correction\n",
+ " cellTable=np.zeros(data.shape[2], np.int32)) \n",
+ " histCalCommonModeCor.fill(data) # filling histogram with common mode corrected data\n",
+ " \n",
+ " if cm_mean_im is None:\n",
+ " cm_mean_im = np.nanmean(data, axis=2) \n",
+ " cm_single_im = data[...,0] # The common mode corrected image corresponding to the first frame \n",
+ " \n",
+ " ddsetcm[...] = np.moveaxis(data, 2, 0)\n",
"\n",
- " patternClassifierLH._noisemap = noise[:, :x//2, :]\n",
- " patternClassifierUH._noisemap = noise[:, x//2:, :]\n",
- " data = cmCorrection.correct(data.astype(np.float32),\n",
- " cellTable=np.zeros(data.shape[2], np.int32))\n",
- " ddsetcm[...] = np.moveaxis(data, 2, 0)\n",
- "\n",
- " dataLH = data[:, :x//2, :]\n",
- " dataUH = data[:, x//2:, :]\n",
+ " # Dividing the data into left and right hemispheres:\n",
+ " dataLH = data[:, :pixels_x//2, :]\n",
+ " dataRH = data[:, pixels_x//2:, :]\n",
"\n",
- " dataLH, patternsLH = patternClassifierLH.classify(dataLH)\n",
- " dataUH, patternsUH = patternClassifierUH.classify(dataUH)\n",
+ " dataLH, patternsLH = patternClassifierLH.classify(dataLH) # pattern classification on corrected data\n",
+ " dataRH, patternsRH = patternClassifierRH.classify(dataRH)\n",
"\n",
- " data[:, :x//2, :] = dataLH\n",
- " data[:, x//2:, :] = dataUH\n",
+ " data[:, :pixels_x//2, :] = dataLH\n",
+ " data[:, pixels_x//2:, :] = dataRH\n",
"\n",
" patterns = np.zeros(data.shape, patternsLH.dtype)\n",
- " patterns[:, :x//2, :] = patternsLH\n",
- " patterns[:, x//2:, :] = patternsUH\n",
+ " patterns[:, :pixels_x//2, :] = patternsLH\n",
+ " patterns[:, pixels_x//2:, :] = patternsRH\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",
+ " histCalPcorr.fill(data) # filling histogram with split events corrected data\n",
+ " \n",
+ " singles_events = np.zeros(data.shape, data.dtype)\n",
+ " for s in range(100, 102): # events' patterns indices are: 100 (singles), 101 (first singles)\n",
+ " singles_events = copy.copy(data[...])\n",
+ " singles_events[patterns != s] = np.nan # Throwing out doubles, triples, quadruple, clusters\n",
+ " histCalPcorrS.fill(singles_events) # filling histogram with singles events data\n",
+ " del singles_events # This line and the next delet this variable from the memory so we do not \n",
+ " gc.collect() # fill the memory up with copies of data per sequence\n",
+ " \n",
+ " # mean_im_cc = averaged over all pattern classified corrected images \n",
+ " # (in the first sequence only)\n",
" if mean_im_cc is None:\n",
- " mean_im_cc = np.nanmean(data, axis=2)\n",
- " single_im_cc = data[...,0]\n",
- " \n",
+ " mean_im_cc = np.nanmean(data, axis=2) \n",
+ " single_im_cc = data[...,0] # The final corrected image corresponding to the first frame\n",
" except Exception as e:\n",
- " print(\"Couldn't calibrate data in {}: {}\".format(f, e))\n"
+ " print(f\"Couldn't calibrate data in {f}: {e}\\n\")\n",
+ "\n",
+ "print(\"In addition to offset correction, the following corrections were performed:\")\n",
+ "for k, v in corr_bools.items():\n",
+ " if v:\n",
+ " print(k)"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2018-12-06T16:10:56.126409Z",
- "start_time": "2018-12-06T16:10:56.096242Z"
- }
- },
+ "metadata": {},
"outputs": [],
"source": [
- "if do_pattern_classification:\n",
- " print(\"******************LOWER HEMISPHERE******************\\n\")\n",
+ "# Histograming the resulting spectra:\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",
+ "# First _ refers to the bin edges and second _ refers to statistics and we ignore both. \n",
"\n",
- " ax = fig.add_subplot(4,4,2+i)\n",
+ "# if you use histCalRaw.get(cumulatative = True) and so on, the cumulatative = True turns the counts array such as \n",
+ "# RawHistVals and so on into a 1D array instead of keeping the original shape:\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",
+ "RawHistVals, _, RawHistMids, _ = histCalRaw.get() \n",
+ "off_cor_HistVals, _, off_cor_HistMids, _ = histCalOffsetCor.get()\n",
+ "if corr_bools.get('common_mode'):\n",
+ " cm_cor_HistVals, _, cm_HistMids, _ = histCalCommonModeCor.get()\n",
+ "if corr_bools.get('relgain'):\n",
+ " gain_cor_HistVals, _, gain_cor_HistMids, _ = histCalGainCor.get()\n",
+ "split_HistVals, _, split_HistMids, _ = histCalPcorr.get() # split events corrected\n",
+ "singles_HistVals, _, singles_HistMids, _ = histCalPcorrS.get() # last s in variable names: singles events"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Saving intermediate data to disk:\n",
+ "\n",
+ "np.savez(os.path.join(out_folder, 'Raw_Events.npz'), RawHistMids, RawHistVals)\n",
+ "np.savez(os.path.join(out_folder, 'Offset_Corrected_Events.npz'), off_cor_HistMids, off_cor_HistVals)\n",
+ "if corr_bools.get('common_mode'):\n",
+ " np.savez(os.path.join(out_folder, 'Common_Mode_Corrected_Events.npz'), cm_HistMids, cm_cor_HistVals)\n",
+ "if corr_bools.get('relgain'):\n",
+ " np.savez(os.path.join(out_folder, 'Gain_Corrected_Events.npz'), gain_cor_HistMids, gain_cor_HistVals)\n",
+ "np.savez(os.path.join(out_folder, 'Split_Events.npz'), split_HistMids, split_HistVals)\n",
+ "np.savez(os.path.join(out_folder, 'Singles_Events.npz'), singles_HistMids, singles_HistVals)\n",
+ "\n",
+ "print(\"Various spectra are saved to disk in the form of histograms. Please check {}\".format(out_folder))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# depending on the gain, the ranges of data are different, so we set a few parameters so that the plots always look\n",
+ "# good.\n",
"\n",
- " mmaps = [\"doublemap\", \"triplemap\", \"quadmap\"]\n",
- " k = 0\n",
- " for i, m in enumerate(mmaps):\n",
+ "if gain == 1:\n",
+ " x_range = (3000, 16000)\n",
+ "elif gain == 64:\n",
+ " x_range = (0, 1000)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "display(Markdown('### Raw vs. Corrected Spectra'))\n",
+ "\n",
+ "figure = [{'x': RawHistMids,\n",
+ " 'y': RawHistVals,\n",
+ " 'y_err': np.sqrt(RawHistVals[:]),\n",
+ " 'drawstyle': 'steps-post',\n",
+ " 'errorstyle': 'bars',\n",
+ " 'errorcoarsing': 2,\n",
+ " 'label': 'Uncorrected'\n",
+ " },\n",
+ " {'x': off_cor_HistMids,\n",
+ " 'y': off_cor_HistVals,\n",
+ " 'y_err': np.sqrt(off_cor_HistVals[:]),\n",
+ " 'drawstyle': 'steps-post',\n",
+ " 'errorstyle': 'bars',\n",
+ " 'errorcoarsing': 2,\n",
+ " 'label': 'Offset Corrected'\n",
+ " }]\n",
+ "\n",
+ "if corr_bools.get('common_mode'):\n",
+ " figure.append({'x': cm_HistMids,\n",
+ " 'y': cm_cor_HistVals,\n",
+ " 'y_err': np.sqrt(cm_cor_HistVals[:]),\n",
+ " 'drawstyle': 'steps-post',\n",
+ " 'errorstyle': 'bars',\n",
+ " 'errorcoarsing': 2,\n",
+ " 'label': 'Common Mode Corrected'})\n",
+ "\n",
+ "if corr_bools.get('relgain'):\n",
+ " xrange = (cti_limit_low, cti_limit_high)\n",
+ " figure.append({'x': gain_cor_HistMids,\n",
+ " 'y': gain_cor_HistVals,\n",
+ " 'y_err': np.sqrt(gain_cor_HistVals[:]),\n",
+ " 'drawstyle': 'steps-post',\n",
+ " 'errorstyle': 'bars',\n",
+ " 'errorcoarsing': 2,\n",
+ " 'label': 'Gain Corrected'})\n",
+ "else:\n",
+ " xrange = x_range\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"
+ " \n",
+ "if corr_bools.get('pattern_class'): \n",
+ " figure.extend([{'x': split_HistMids,\n",
+ " 'y': split_HistVals,\n",
+ " 'y_err': np.sqrt(split_HistVals[:]),\n",
+ " 'drawstyle': 'steps-post',\n",
+ " 'errorstyle': 'bars',\n",
+ " 'errorcoarsing': 2,\n",
+ " 'label': 'Split Events Corrected'\n",
+ " },\n",
+ " {'x': singles_HistMids,\n",
+ " 'y': singles_HistVals,\n",
+ " 'y_err': np.sqrt(singles_HistVals[:]),\n",
+ " 'drawstyle': 'steps-post',\n",
+ " 'errorstyle': 'bars',\n",
+ " 'errorcoarsing': 2,\n",
+ " 'label': 'Singles Events'\n",
+ " }])\n",
+ "fig = xana.simplePlot(figure, aspect=1, x_label='ADU', y_label='Number of Occurrences', figsize='2col',\n",
+ " y_log=True, x_range=xrange, title = '1 ADU per bin is used.',\n",
+ " legend='top-right-frame-1col')"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2018-12-06T16:10:56.176160Z",
- "start_time": "2018-12-06T16:10:56.127853Z"
- }
- },
+ "metadata": {},
"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",
+ "# This function plots pattern statistics:\n",
+ "\n",
+ "def classification_plot(patternStats, hemisphere):\n",
+ " \n",
+ " print(\"****************** {} HEMISPHERE ******************\\n\"\n",
+ " .format(hemisphere))\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",
+ " relativeOccurances = []\n",
" labels = []\n",
" for i, f in enumerate(sfields):\n",
- " relativeOccurances.append(patternStatsUH[f])\n",
+ " relativeOccurances.append(patternStats[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",
+ " for k in range(len(patternStats[f])):\n",
+ " relativeOccurances.append(patternStats[f][k])\n",
+ " labels.append(\"{}({})\".format(f, k))\n",
" relativeOccurances = np.array(relativeOccurances, np.float)\n",
- " relativeOccurances/=np.sum(relativeOccurances)\n",
+ " relativeOccurances /= np.sum(relativeOccurances)\n",
" pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)\n",
- " ax.set_title(\"Pattern occurrence\")\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 = fig.add_subplot(4, 4, 2+i)\n",
+ " pmap = ax.imshow(patternStats[m], interpolation=\"nearest\", vmax=2*np.nanmedian(patternStats[m]))\n",
" ax.set_title(m)\n",
" cb = fig.colorbar(pmap)\n",
"\n",
@@ -786,44 +934,64 @@
" 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",
+ " ax = fig.add_subplot(4, 4, 2+len(smaps)+k)\n",
+ " pmap = ax.imshow(patternStats[m][j], interpolation=\"nearest\", vmax=2*np.median(patternStats[m][j]))\n",
+ " ax.set_title(\"{}({})\".format(m,j))\n",
" cb = fig.colorbar(pmap)\n",
- " k+=1\n",
- "\n",
- " print(\"******************UPPER HEMISPHERE******************\\n\") "
+ " k+=1"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2018-12-06T16:10:56.190150Z",
- "start_time": "2018-12-06T16:10:56.177570Z"
- }
- },
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# The next two cells plot the classification results for left and right hemispheres, respectively:\n",
+ "display(Markdown('### Classification Results - Plots'))\n",
+ "if corr_bools.get('pattern_class'):\n",
+ " patternStatsLH = patternClassifierLH.getPatternStats()\n",
+ " classification_plot(patternStatsLH, 'Left')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
"outputs": [],
"source": [
- "if do_pattern_classification:\n",
+ "if corr_bools.get('pattern_class'):\n",
+ " patternStatsRH = patternClassifierRH.getPatternStats()\n",
+ " classification_plot(patternStatsRH, 'Right')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "display(Markdown('### Classification Results - Tabulated Statistics'))\n",
+ "\n",
+ "if corr_bools.get('pattern_class'):\n",
" t0 = PrettyTable()\n",
- " t0.title = \"Total number of Counts after all corrections\"\n",
- " t0.field_names = [\"Hemisphere\",\"Singles\", \"First-singles\", \"Clusters\"]\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",
+ " t0.add_row([\"RH\", patternStatsRH['singles'], patternStatsRH['first singles'], patternStatsRH['clusters']])\n",
" print(t0)\n",
- "\n",
+ " \n",
+ " print(\"Abbreviations: D (Doubles), T (Triples), Q (Quadruples), L (Left), R (Right), and H (Hemisphere).\")\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",
+ " t1.field_names = [\"Index\", \"D-LH\", \"D-RH\", \"T-LH\", \"T-RH\", \"Q-LH\", \"Q-RH\"]\n",
+ " t1.add_row([0, patternStatsLH['doubles'][0], patternStatsRH['doubles'][0], patternStatsLH['triples'][0], \n",
+ " patternStatsRH['triples'][0], patternStatsLH['quads'][0], patternStatsRH['quads'][0]])\n",
+ " t1.add_row([1, patternStatsLH['doubles'][1], patternStatsRH['doubles'][1], patternStatsLH['triples'][1], \n",
+ " patternStatsRH['triples'][1], patternStatsLH['quads'][1], patternStatsRH['quads'][1]])\n",
+ " t1.add_row([2, patternStatsLH['doubles'][2], patternStatsRH['doubles'][2], patternStatsLH['triples'][2], \n",
+ " patternStatsRH['triples'][2], patternStatsLH['quads'][2], patternStatsRH['quads'][2]])\n",
+ " t1.add_row([3, patternStatsLH['doubles'][3], patternStatsRH['doubles'][3], patternStatsLH['triples'][3], \n",
+ " patternStatsRH['triples'][3], patternStatsLH['quads'][3], patternStatsRH['quads'][3]])\n",
" print(t1)"
]
},
@@ -832,33 +1000,39 @@
"execution_count": null,
"metadata": {
"ExecuteTime": {
- "end_time": "2018-12-06T16:10:56.203219Z",
- "start_time": "2018-12-06T16:10:56.191509Z"
+ "end_time": "2018-12-06T16:10:56.190150Z",
+ "start_time": "2018-12-06T16:10:56.177570Z"
}
},
"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",
+ "if corr_bools.get('pattern_class'):\n",
+ " doublesLH = patternStatsLH['doubles'][0] + patternStatsLH['doubles'][1] + patternStatsLH['doubles'][2] + \\\n",
+ " patternStatsLH['doubles'][3]\n",
+ " triplesLH = patternStatsLH['triples'][0] + patternStatsLH['triples'][1] + patternStatsLH['triples'][2] + \\\n",
+ " patternStatsLH['triples'][3]\n",
+ " quadsLH = patternStatsLH['quads'][0] + patternStatsLH['quads'][1] + patternStatsLH['quads'][2] + \\\n",
+ " 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",
+ " doublesRH = patternStatsRH['doubles'][0] + patternStatsRH['doubles'][1] + patternStatsRH['doubles'][2] + \\\n",
+ " patternStatsRH['doubles'][3]\n",
+ " triplesRH = patternStatsRH['triples'][0] + patternStatsRH['triples'][1] + patternStatsRH['triples'][2] + \\\n",
+ " patternStatsRH['triples'][3]\n",
+ " quadsRH = patternStatsRH['quads'][0] + patternStatsRH['quads'][1] + patternStatsRH['quads'][2] + \\\n",
+ " patternStatsRH['quads'][3]\n",
+ " allsinglesRH = patternStatsRH['singles'] + patternStatsRH['first singles']\n",
+ " eventsRH = allsinglesRH + doublesRH + triplesRH + quadsRH\n",
"\n",
" if eventsLH > 0.:\n",
" reloccurLH = np.array([allsinglesLH/eventsLH, doublesLH/eventsLH, triplesLH/eventsLH, quadsLH/eventsLH])\n",
" else:\n",
" reloccurLH = np.array([0]*4)\n",
- " if eventsUH > 0.:\n",
- " reloccurUH = np.array([allsinglesUH/eventsUH, doublesUH/eventsUH, triplesUH/eventsUH, quadsUH/eventsUH])\n",
+ " if eventsRH > 0.:\n",
+ " reloccurRH = np.array([allsinglesRH/eventsRH, doublesRH/eventsRH, triplesRH/eventsRH, quadsRH/eventsRH])\n",
" else:\n",
- " reloccurUH = np.array([0]*4)"
+ " reloccurRH = np.array([0]*4)"
]
},
{
@@ -866,32 +1040,63 @@
"execution_count": null,
"metadata": {
"ExecuteTime": {
- "end_time": "2018-12-06T16:10:56.212586Z",
- "start_time": "2018-12-06T16:10:56.204731Z"
+ "end_time": "2018-12-06T16:10:56.203219Z",
+ "start_time": "2018-12-06T16:10:56.191509Z"
}
},
"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",
+ "display(Markdown('### Classification Results - Pie Charts'))\n",
+ "\n",
+ "if corr_bools.get('pattern_class'):\n",
+ " fig = plt.figure(figsize=(12, 7))\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",
+ " ax.set_title(\"Pattern Occurrence in 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",
+ " ax = fig.add_subplot(1, 2, 2)\n",
+ " pie = ax.pie(reloccurRH, labels=labels, autopct='%1.1f%%', shadow=True)\n",
+ " ax.set_title(\"Pattern Occurrence in RH\")\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:10:56.212586Z",
+ "start_time": "2018-12-06T16:10:56.204731Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "# Again, we are setting a few parameters so that the plots always look good regardless of what the gain is:\n",
+ "\n",
+ "if gain == 1:\n",
+ " vmin = -50 \n",
+ " vmax = 40000\n",
+ " panel_top_low_lim = -3000\n",
+ " panel_top_high_lim = 65000\n",
+ " panel_side_low_lim = -3000 \n",
+ " panel_side_high_lim = 65000\n",
+ "elif gain == 64:\n",
+ " vmin = 0\n",
+ " vmax = 20000\n",
+ " panel_top_low_lim = 0\n",
+ " panel_top_high_lim = 20000\n",
+ " panel_side_low_lim = 0\n",
+ " panel_side_high_lim = 20000"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Mean Image of first Sequence ##"
+ "### Mean Images of the First Sequence ###"
]
},
{
@@ -899,33 +1104,77 @@
"execution_count": null,
"metadata": {
"ExecuteTime": {
- "end_time": "2018-12-06T16:11:08.317130Z",
- "start_time": "2018-12-06T16:11:05.788655Z"
- }
+ "end_time": "2018-12-06T16:10:56.212586Z",
+ "start_time": "2018-12-06T16:10:56.204731Z"
+ },
+ "scrolled": false
},
"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=70000)\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=70000)"
+ "fig = xana.heatmapPlot(uncor_mean_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1, \n",
+ " x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=vmin, vmax=vmax, \n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " panel_top_low_lim = panel_top_low_lim, panel_top_high_lim = panel_top_high_lim,\n",
+ " panel_side_low_lim = panel_side_low_lim, panel_side_high_lim = panel_side_high_lim, \n",
+ " title = 'Uncorrected Image Averaged over the First Sequence')\n",
+ "\n",
+ "# At lower gain, we need to change the vmin and vmax for the following plots:\n",
+ "if gain == 64:\n",
+ " vmin = 0\n",
+ " vmax = 8000\n",
+ " \n",
+ "\n",
+ "fig = xana.heatmapPlot(offset_mean_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1, \n",
+ " x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=vmin, vmax=vmax, \n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " panel_top_low_lim = panel_top_low_lim, panel_top_high_lim = panel_top_high_lim,\n",
+ " panel_side_low_lim = panel_side_low_lim, panel_side_high_lim = panel_side_high_lim, \n",
+ " title = 'Offset Corrected Image Averaged over the First Sequence')\n",
+ "if corr_bools.get('relgain'):\n",
+ " fig = xana.heatmapPlot(gain_mean_im, x_label='Columns', y_label='Rows', \n",
+ " lut_label='Gain Corrected Signal (ADU)', aspect=1, \n",
+ " x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=vmin, vmax=vmax, \n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " panel_top_low_lim = panel_top_low_lim, panel_top_high_lim = panel_top_high_lim,\n",
+ " panel_side_low_lim = panel_side_low_lim, panel_side_high_lim = panel_side_high_lim, \n",
+ " title = 'Mean Gain Corrected Image of the First Sequence')\n",
+ "\n",
+ "if corr_bools.get('pattern_class') and corr_bools.get('common_mode'):\n",
+ " fig = xana.heatmapPlot(cm_mean_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1, \n",
+ " x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=vmin, vmax=vmax,\n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " panel_top_low_lim = panel_top_low_lim, panel_top_high_lim = panel_top_high_lim,\n",
+ " panel_side_low_lim = panel_side_low_lim, panel_side_high_lim = panel_side_high_lim, \n",
+ " title = 'Common Mode Corrected Image Averaged over the First Sequence')\n",
+ " \n",
+ "# Again, we need to change the vmin and vmax and other settings for this last plot:\n",
+ " if gain == 1:\n",
+ " vmin = 0\n",
+ " vmax = 600\n",
+ " panel_top_low_lim = 0\n",
+ " panel_top_high_lim = 500\n",
+ " panel_side_low_lim = 0\n",
+ " panel_side_high_lim = 2000\n",
+ " elif gain == 64:\n",
+ " vmin = 0\n",
+ " vmax = 600\n",
+ " panel_top_low_lim = 0\n",
+ " panel_top_high_lim = 2000\n",
+ " panel_side_low_lim = 0\n",
+ " panel_side_high_lim = 2000\n",
+ " \n",
+ " fig = xana.heatmapPlot(mean_im_cc, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,\n",
+ " x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=vmin, vmax=vmax,\n",
+ " panel_top_low_lim = panel_top_low_lim, panel_top_high_lim = panel_top_high_lim,\n",
+ " panel_side_low_lim = panel_side_low_lim, panel_side_high_lim = panel_side_high_lim, \n",
+ " title = 'Image after All Corrections Averaged over the First Sequence')"
]
},
{
"cell_type": "markdown",
- "metadata": {
- "collapsed": true
- },
+ "metadata": {},
"source": [
- "## Single Shot of first Sequnce ##"
+ "### Images of the First Frame of the First Sequence ###"
]
},
{
@@ -933,104 +1182,196 @@
"execution_count": null,
"metadata": {
"ExecuteTime": {
- "end_time": "2018-12-06T16:11:10.908912Z",
- "start_time": "2018-12-06T16:11:08.318486Z"
- }
+ "end_time": "2018-12-06T16:11:08.317130Z",
+ "start_time": "2018-12-06T16:11:05.788655Z"
+ },
+ "scrolled": false
},
"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=70000)\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=70000)"
+ "if gain == 1:\n",
+ " vmin = -50\n",
+ " vmax = 10000\n",
+ " panel_top_low_lim = 0\n",
+ " panel_top_high_lim = 20000\n",
+ " panel_side_low_lim = 0\n",
+ " panel_side_high_lim = 20000\n",
+ " \n",
+ "elif gain == 64:\n",
+ " vmin = 0\n",
+ " vmax = 2000\n",
+ " panel_top_low_lim = 0\n",
+ " panel_top_high_lim = 20000\n",
+ " panel_side_low_lim = 0\n",
+ " panel_side_high_lim = 20000\n",
+ "\n",
+ "fig = xana.heatmapPlot(uncor_single_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1, \n",
+ " x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=vmin, vmax=vmax,\n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " panel_top_low_lim = panel_top_low_lim, panel_top_high_lim = panel_top_high_lim,\n",
+ " panel_side_low_lim = panel_side_low_lim, panel_side_high_lim = panel_side_high_lim, \n",
+ " title = 'Uncorrected Image (First Frame of First Sequence)')\n",
+ "\n",
+ "fig = xana.heatmapPlot(offset_single_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1, \n",
+ " x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=vmin, vmax=vmax,\n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " panel_top_low_lim = panel_top_low_lim, panel_top_high_lim = panel_top_high_lim,\n",
+ " panel_side_low_lim = panel_side_low_lim, panel_side_high_lim = panel_side_high_lim, \n",
+ " title = 'Offset Corrected Image (First Frame of First Sequence)')\n",
+ "\n",
+ "if corr_bools.get('pattern_class'):\n",
+ " # At lower gain, we need to change the vmin and vmax for the following plots:\n",
+ " if gain == 64:\n",
+ " vmin = 0\n",
+ " vmax = 600\n",
+ " panel_top_low_lim = 0\n",
+ " panel_top_high_lim = 800\n",
+ " panel_side_low_lim = 0\n",
+ " panel_side_high_lim = 800\n",
+ " if corr_bools.get('common_mode'):\n",
+ " fig = xana.heatmapPlot(cm_single_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1, \n",
+ " x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=vmin, vmax=vmax,\n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " panel_top_low_lim = panel_top_low_lim, panel_top_high_lim = panel_top_high_lim,\n",
+ " panel_side_low_lim = panel_side_low_lim, panel_side_high_lim = panel_side_high_lim, \n",
+ " title = 'Common Mode Corrected Image (First Frame of First Sequence)')\n",
+ " \n",
+ " fig = xana.heatmapPlot(single_im_cc, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,\n",
+ " x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=vmin, vmax=vmax,\n",
+ " panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)', \n",
+ " panel_top_low_lim = panel_top_low_lim, panel_top_high_lim = panel_top_high_lim,\n",
+ " panel_side_low_lim = panel_side_low_lim, panel_side_high_lim = panel_side_high_lim, \n",
+ " title = 'Image after All Corrections (First Frame of First Sequence)')"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "scrolled": false
- },
+ "metadata": {},
"outputs": [],
"source": [
- "single = None\n",
- "doubles = []\n",
- "triples = []\n",
- "quads = []\n",
- "import copy\n",
- "with h5py.File(\"{}/CORR-R{:04d}-PNCCD01-S{:05d}.h5\".format(out_folder, run, sequences[0]), 'r', driver='core') as infile:\n",
- " data = infile[h5path+\"/pixels_classified\"][()].astype(np.float)\n",
- " patterns = infile[h5path+\"/patterns\"][()].astype(np.float)\n",
- "\n",
- " i = 4\n",
- " single = copy.copy(data[...])\n",
- " single[patterns != 100] = np.nan\n",
- " for d in range(200,204):\n",
- " double = copy.copy(data[...])\n",
- " double[patterns != d] = np.nan\n",
- " doubles.append(double)\n",
- " for d in range(300,304):\n",
- " triple = copy.copy(data[...])\n",
- " triple[patterns != d] = np.nan\n",
- " triples.append(triple)\n",
- " \n",
- " for d in range(400,404):\n",
- " quad = copy.copy(data[...])\n",
- " quad[patterns != d] = np.nan\n",
- " quads.append(quad)"
+ "# Resetting the histogram calculators:\n",
+ "histCalRaw.reset()\n",
+ "histCalOffsetCor.reset()\n",
+ "histCalPcorr.reset()\n",
+ "histCalPcorrS.reset()\n",
+ "if corr_bools.get('common_mode'):\n",
+ " histCalCommonModeCor.reset()\n",
+ "if corr_bools.get('relgain'):\n",
+ " histCalGainCor.reset()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next, the corrected event patterns are read from the patterns/ dataset created previously and are separated into 4 different categories (singles, doubles, triples and quadruples) using the pattern indices. However, this is done only for one sequence, corresponding to the seq_num variable, as an example. \n",
+ "\n",
+ "Note that the number of bins and the bin range for the following histograms may be different from those presented above (depending on gain) to make the counts more noticible and the peaks more defined.\n",
+ "\n",
+ "If you are interested in plotting the events from all sequences or the spectra of half of the sensor, execute the spectra_pnCCD_NBC.ipynb notebook."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "scrolled": false
- },
+ "metadata": {},
"outputs": [],
"source": [
- "hA = 0\n",
- "h, e = np.histogram(single.flatten(), bins=1000, range=(-500, 50000))\n",
- "fig = plt.figure(figsize=(10,7))\n",
- "ax = fig.add_subplot(111)\n",
- "ax.step(e[:-1], h, color='blue', label='single pixel events')\n",
- "ax.semilogy()\n",
- "ax.set_xlabel(\"Energy (ADU)\")\n",
- "ax.set_ylabel(\"Events\")\n",
- "hA += h\n",
- "h = 0\n",
- "for double in doubles:\n",
- " hd, e = np.histogram(double.flatten(), bins=1000, range=(-500, 50000))\n",
- " h += hd\n",
- " hA += hd\n",
- " \n",
- "ax.step(e[:-1], h, color='red', label='double split events')\n",
+ "if corr_bools.get('pattern_class'):\n",
+ " singles = []\n",
+ " doubles = []\n",
+ " triples = []\n",
+ " quads = []\n",
"\n",
- "h = 0\n",
- "for triple in triples:\n",
- " hd, e = np.histogram(triple.flatten(), bins=1000, range=(-500, 50000))\n",
- " h += hd\n",
- " hA += hd\n",
- " \n",
- "ax.step(e[:-1], h, color='green', label='triple split events')\n",
+ " with h5py.File(\"{}/CORR-R{:04d}-PNCCD01-S{:05d}.h5\".format(out_folder, run, sequences[seq_num]), 'r', \n",
+ " driver='core') as infile:\n",
"\n",
- "h = 0\n",
- "for quad in quads:\n",
- " hd, e = np.histogram(quad.flatten(), bins=1000, range=(-500, 50000))\n",
- " h += hd\n",
- " hA += hd\n",
- " \n",
- "ax.step(e[:-1], h, color='purple', label='quadruple split events')\n",
- "ax.step(e[:-1], hA, color='grey', label='all events')\n",
+ " data = infile[h5path+\"/pixels_classified\"][()].astype(np.float32) # classifications\n",
+ " patterns = infile[h5path+\"/patterns\"][()].astype(np.float32) # event patterns\n",
+ " \n",
+ " # events' patterns indices are as follows: 100 (singles), 101 (first singles), 200 - 203 (doubles),\n",
+ " # 300 - 303 (triples), and 400 - 403 (quadruples). Note that for the last three types of patterns, \n",
+ " # there are left, right, up, and down indices.\n",
+ "\n",
+ " # Separating the events:\n",
+ " # Singles and First Singles:\n",
+ " for s in range(100, 102):\n",
+ " single = copy.copy(data[...])\n",
+ " single[patterns != s] = np.nan\n",
+ " singles.append(single)\n",
+ " \n",
+ "\n",
+ " for d in range(200, 204):\n",
+ " double = copy.copy(data[...])\n",
+ " double[patterns != d] = np.nan\n",
+ " doubles.append(double)\n",
"\n",
- "l = ax.legend()"
+ " for t in range(300, 304):\n",
+ " triple = copy.copy(data[...])\n",
+ " triple[patterns != t] = np.nan\n",
+ " triples.append(triple) \n",
+ "\n",
+ " for q in range(400, 404):\n",
+ " quad = copy.copy(data[...])\n",
+ " quad[patterns != q] = np.nan\n",
+ " quads.append(quad)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "if corr_bools.get('pattern_class'):\n",
+ " hA = 0\n",
+ " h = 0\n",
+ " for single in singles:\n",
+ " hs, e = np.histogram(single.flatten(), bins=event_bins, range=b_range) # h: histogram counts, e: bin edges\n",
+ " h += hs\n",
+ " hA += hs # hA: counts all events (see below)\n",
+ "\n",
+ " \n",
+ " # bin edges array has one extra element => need to plot from 0 to the one before the last element to have the \n",
+ " # same size as h-array => in what follows, we use e[:-1] (-1 means one before the last element)\n",
+ "\n",
+ " #TODO adapt the title depending on corrections\n",
+ " display(Markdown('### Histograms of Offset, Common Mode and Gain Corrected Events for One Sequence Only'))\n",
+ " fig = plt.figure(figsize=(10, 7))\n",
+ " ax = fig.add_subplot(111)\n",
+ " ax.step(e[:-1], h, color='blue', label='Events Involving Single Pixels Only') \n",
+ " ax.semilogy() # y-axis is log, x-axis is linear\n",
+ " ax.set_xlabel(\"Energy (ADU) [{} bins per {} ADU.]\".format(event_bins, b_range[1]-b_range[0]))\n",
+ " ax.set_ylabel(\"Gain Corrected Events (counts) for One Sequence\")\n",
+ " ax.set_xlim(x_range)\n",
+ "\n",
+ " h = 0\n",
+ " for double in doubles:\n",
+ " hd, e = np.histogram(double.flatten(), bins=event_bins, range=b_range) \n",
+ " h += hd\n",
+ " hA += hd\n",
+ "\n",
+ " ax.step(e[:-1], h, color='red', label='Events Splitting on Double Pixels')\n",
+ "\n",
+ " h = 0\n",
+ " for triple in triples:\n",
+ " ht, e = np.histogram(triple.flatten(), bins=event_bins, range=b_range) \n",
+ " h += ht\n",
+ " hA += ht\n",
+ "\n",
+ " ax.step(e[:-1], h, color='green', label='Events Splitting on Triple Pixels')\n",
+ "\n",
+ " h = 0\n",
+ " for quad in quads:\n",
+ " hq, e = np.histogram(quad.flatten(), bins=event_bins, range=b_range) \n",
+ " h += hq\n",
+ " hA += hq\n",
+ "\n",
+ " ax.step(e[:-1], h, color='purple', label='Events Splitting on Quadruple Pixels')\n",
+ "\n",
+ " ax.step(e[:-1], hA, color='grey', label='All Valid Events')\n",
+ " l = ax.legend()"
]
},
{
diff --git a/xfel_calibrate/notebooks.py b/xfel_calibrate/notebooks.py
index 03053917b..12c97b96a 100644
--- a/xfel_calibrate/notebooks.py
+++ b/xfel_calibrate/notebooks.py
@@ -96,11 +96,16 @@ notebooks = {
"default concurrency": None,
"cluster cores": 32},
},
+ "RELGAIN": {
+ "notebook": "notebooks/pnCCD/Characterize_pnCCD_Gain.ipynb",
+ "concurrency": {"parameter": None,
+ "default concurrency": None,
+ "cluster cores": 32},
+ },
"CORRECT": {
"notebook": "notebooks/pnCCD/Correct_pnCCD_NBC.ipynb",
- "concurrency": {"parameter": "sequences",
- "default concurrency": [-1],
- "use function": "balance_sequences",
+ "concurrency": {"parameter": None,
+ "default concurrency": None,
"cluster cores": 32},
},
},
--
GitLab