Skip to content
Snippets Groups Projects

[LPD][Dark] fix: Remove hard check for consistent gain pixels

Merged [LPD][Dark] fix: Remove hard check for consistent gain pixels
fix/remove_hard_check into master
All threads resolved!
Merged Karim Ahmed requested to merge fix/remove_hard_check into master
All threads resolved!
Compare
and
Show latest version
1 file
+ 4
− 1
Compare changes
  • Side-by-side
  • Inline
notebooks/LPD/LPDChar_Darks_NBC.ipynb
+ 4
− 1
%% Cell type:markdown id: tags:
# LPD Offset, Noise and Dead Pixels Characterization #
Author: M. Karnevskiy, S. Hauf
This notebook performs re-characterize of dark images to derive offset, noise and bad-pixel maps. All three types of constants are evaluated per-pixel and per-memory cell.
The notebook will correctly handle veto settings, but note that if you veto cells you will not be able to use these offsets for runs with different veto settings - vetoed cells will have zero offset.
The evaluated calibration constants are stored locally and injected in the calibration data base.
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/FXE/202304/p003338/raw" # path to input data, required
out_folder = "/gpfs/exfel/data/scratch/kluyvert/lpd-darks-p3338-r133-134-135/" # path to output to, required
metadata_folder = "" # Directory containing calibration_metadata.yml when run by xfel-calibrate
modules = [-1] # list of modules to evaluate, RANGE ALLOWED
run_high = 133 # run number in which high gain data was recorded
run_med = 134 # run number in which medium gain data was recorded
run_low = 135 # run number in which low gain data was recorded
run = 88
operation_mode = 'LPD_ADAPTIVE_GAIN'
karabo_id = "FXE_DET_LPD1M-1" # karabo karabo_id
karabo_da = ['-1'] # a list of data aggregators names, Default [-1] for selecting all data aggregators
source_name = "{}/DET/{}CH0:xtdf" # Source name for raw detector data - filled with karabo_id & module number
ctrl_src_template = "{}/COMP/FEM_MDL_COMP" # Control device source name template.
use_dir_creation_date = True # use the creation date of the directory for database time derivation
cal_db_interface = "tcp://max-exfl-cal001:8015#8025" # the database interface to use
cal_db_timeout = 300000 # timeout on caldb requests"
local_output = True # output constants locally
db_output = False # output constants to database
capacitor_setting = 5 # capacitor_setting for which data was taken
mem_cells = 512 # number of memory cells used
bias_voltage = 250 # detector bias voltage
thresholds_offset_sigma = 3. # bad pixel relative threshold in terms of n sigma offset
thresholds_offset_hard = [400, 1500] # bad pixel hard threshold
thresholds_noise_sigma = 7. # bad pixel relative threshold in terms of n sigma noise
thresholds_noise_hard = [1, 35] # bad pixel hard threshold
ntrains = 500 # maximum number of trains to use in each gain stage
skip_first_ntrains = 10 # Number of first trains to skip
min_trains = 370 # minimum number of trains needed for each gain stage
drop_last_frames_parallelgain = 1 # Discard last N frames of each gain stage in parallel gain mode
bad_gain_tolerance_parallelgain = 0.2 # Error if more than this proportion of pixels are in the wrong gain stage in parallel gain mode
# Parameters for plotting
skip_plots = False # exit after writing corrected files
high_res_badpix_3d = False # plot bad-pixel summary in high resolution
test_for_normality = False # permorm normality test
inject_cell_order = 'auto' # Include memory cell order as part of the detector condition: auto/always/never
```
%% Cell type:code id: tags:
``` python
import copy
import multiprocessing
import warnings
from datetime import datetime
from functools import partial
from logging import warning
from pathlib import Path
warnings.filterwarnings('ignore')
import dateutil.parser
import matplotlib
import pasha as psh
import scipy.stats
from IPython.display import Latex, Markdown, display
matplotlib.use("agg")
import matplotlib.patches as patches
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import tabulate
import yaml
from iCalibrationDB import Conditions, Constants, Detectors, Versions
from XFELDetAna.plotting.heatmap import heatmapPlot
from XFELDetAna.plotting.simpleplot import simplePlot
from extra_data import RunDirectory, by_id
from cal_tools.enums import BadPixels
from cal_tools.lpdlib import (
make_cell_order_condition,
sort_dark_runs_by_gain,
)
from cal_tools.plotting import (
create_constant_overview,
plot_badpix_3d,
show_overview,
show_processed_modules,
)
from cal_tools.tools import (
get_dir_creation_date,
get_from_db,
get_notebook_name,
get_pdu_from_db,
get_random_db_interface,
get_report,
map_gain_stages,
module_index_to_qm,
parse_runs,
raw_data_location_string,
reorder_axes,
run_prop_seq_from_path,
save_const_to_h5,
send_to_db,
)
```
%% Cell type:code id: tags:
``` python
gains = np.arange(3)
max_cells = mem_cells
cells = np.arange(max_cells)
gain_names = ['High', 'Medium', 'Low']
is_adaptive_gain = False
is_parallel_gain = False
if operation_mode == 'LPD_ADAPTIVE_GAIN':
is_adaptive_gain = True
# Check dark runs order and sort if needed.
run_nums = sort_dark_runs_by_gain(
raw_folder=in_folder,
runs=[run_high, run_med, run_low],
ctrl_src=ctrl_src_template.format(karabo_id))
elif operation_mode == 'LPD_PARALLEL_GAIN':
is_parallel_gain = True
run_nums = [run, run, run]
else:
raise ValueError('Invalid LPD operation mode')
for run_num in set(run_nums): # Avoid visiting identical runs more than once.
dc = RunDirectory(Path(in_folder, f'r{run_num:04d}'))
try:
gain_override = bool(dc[ctrl_src_template.format(karabo_id)].run_value('femAsicGainOverride'))
except KeyError:
# Data taken before August 2024 is missing femAsicGainOverride on the FEM
# control device, prevent using parallel gain mode on this data. There are
# only a very small number of parallel gain test runs from before that date.
if is_parallel_gain:
raise ValueError('Control device is lacking data to verify detector '
'is operating in parallel gain mode') from None
continue
if gain_override != is_parallel_gain:
gain_str = 'parallel gain' if gain_override else 'adaptive gain'
raise ValueError(f'Mismatch between chosen operation mode {operation_mode} '
f'and control data indicating {gain_str}') from None
if karabo_da[0] == '-1':
if modules[0] == -1:
modules = list(range(16))
karabo_da = ['LPD{:02d}'.format(i) for i in modules]
else:
modules = [int(x[-2:]) for x in karabo_da]
capacitor_setting_s = f'{capacitor_setting}pf'
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run_nums[0])
print("Using {} as creation time".format(creation_time))
if inject_cell_order not in {'auto', 'always', 'never'}:
raise ValueError("inject_cell_order must be auto/always/never")
run, prop, seq = run_prop_seq_from_path(in_folder)
cal_db_interface = get_random_db_interface(cal_db_interface)
display(Markdown('## Evaluated parameters'))
print('CalDB Interface {}'.format(cal_db_interface))
print("Proposal: {}".format(prop))
print("Memory cells: {}/{}".format(mem_cells, max_cells))
print("Run(s): {}".format(run_nums))
print("Using DB: {}".format(db_output))
print("Input: {}".format(in_folder))
print("Output: {}".format(out_folder))
print("Bias voltage: {}V".format(bias_voltage))
print(f"Capacitor setting: {capacitor_setting_s}")
print(f'Parallel gain: {is_parallel_gain}')
```
%% Cell type:markdown id: tags:
## Data processing
%% Cell type:code id: tags:
``` python
def find_common_cell_pattern(srcdata):
"""The memory cell pattern is not always entirely consistent (seen in p6936 r52, parallel gain mode)
Identify the most common cell pattern, and select data matching that.
"""
cell_ids = srcdata['image.cellId'].drop_empty_trains()[skip_first_ntrains:]
tids_by_cell_pattern = {}
# Count up cell ID patterns per train
for tid, cid_arr in cell_ids.trains():
if cid_arr.ndim > 1:
cid_arr = cid_arr[:, 0]
if is_parallel_gain:
if len(cid_arr) % 3 != 0:
raise ValueError(f"{len(cid_arr)} frames in train {tid} is not divisible by 3 "
"(expected for parallel gain mode)")
# Only the high-gain cell IDs are valid
cid_arr = cid_arr[:len(cid_arr) // 3]
if drop_last_frames_parallelgain:
cid_arr = cid_arr[:-drop_last_frames_parallelgain]
tids_by_cell_pattern.setdefault(tuple(cid_arr), []).append(tid)
if not tids_by_cell_pattern:
raise ValueError("No cell ID patterns were found.")
warning("No cell ID patterns were found.")
return None, None
most_common, sel_tids = max(tids_by_cell_pattern.items(), key=lambda p: len(p[1]))
if len(sel_tids) <= (0.5 * len(cell_ids.train_ids)):
raise ValueError("No cell ID pattern matched over half of the trains "
f"(max {len(sel_tids)} / {len(cell_ids.train_ids)}")
return np.array(most_common), sel_tids[:ntrains]
```
%% Cell type:code id: tags:
``` python
parallel_num_procs = min(6, len(modules)*3)
parallel_num_threads = multiprocessing.cpu_count() // parallel_num_procs
# the actual characterization
def characterize_module(run_path, channel, gg):
run = RunDirectory(run_path, parallelize=False)
det_source = run[source_name.format(karabo_id, channel)]
cellid_pattern, train_ids = find_common_cell_pattern(det_source)
if cellid_pattern is None:
return None, None, None, None, None, None, None, None
data = det_source['image.data'][by_id[train_ids]]
cell_ids = det_source['image.cellId'][by_id[train_ids]]
# If there is no data available, return and expect this
# module to be skipped later.
if len(data.train_ids) == 0:
return None, None, None, None, None, None, None, None
elif len(data.train_ids) < min_trains:
raise Exception(f"Run {run_path} only contains {len(data.train_ids)} trains, but {min_trains} required")
im = data.ndarray()
if im.ndim > 3:
im = im[:, 0] # Drop extra dimension
cellid = cell_ids.ndarray()
if cellid.ndim > 1:
cellid = cellid[:, 0]
# Split off gain bits from pixel data.
gains = (im & 0x3000)
gains >>= 12
im &= 0b0000111111111111
im = im.astype(np.float32)
if is_parallel_gain:
# (trains, gains, frames)
frame_indexes = np.arange(len(cellid)).reshape(len(train_ids), 3, -1)
# The last frame in each gain stage is often bad, so we may drop it
if drop_last_frames_parallelgain:
frame_indexes = frame_indexes[..., :-drop_last_frames_parallelgain]
# Cell IDs from high gain, data from the current gain stage
cellid = cellid[frame_indexes[:, 0].ravel()]
gains = gains[frame_indexes[:, gg].ravel()]
im = im[frame_indexes[:, gg].ravel()]
# Discard frames with too many pixels in the wrong gain stage,
# and check we still have data for each cell ID.
INVALID_CELLID = 9999
wrong_gain_fraction = (gains != gg).mean(axis=(1, 2)) # fraction per frame
cellid[wrong_gain_fraction > bad_gain_tolerance_parallelgain] = INVALID_CELLID
cellids_missing = cellid_pattern[~np.isin(cellid_pattern, cellid)]
if len(cellids_missing):
raise ValueError(f"No frames left for cells {cellids_missing} "
"after discarding frames with wrong gain stage")
im = reorder_axes(im,
from_order=('frames', 'slow_scan', 'fast_scan'),
to_order=('fast_scan', 'slow_scan', 'frames'),
)
context = psh.context.ThreadContext(num_workers=parallel_num_threads)
offset = context.alloc(shape=(im.shape[0], im.shape[1], max_cells), dtype=np.float64)
noise = context.alloc(like=offset)
normal_test = context.alloc(like=offset)
def process_cell(worker_id, array_index, cc):
idx = cellid == cc
im_slice = im[..., idx]
if np.any(idx):
offset[..., cc] = np.median(im_slice, axis=2)
noise[..., cc] = np.std(im_slice, axis=2)
if test_for_normality:
_, normal_test[..., cc] = scipy.stats.normaltest(
im[:, :, idx], axis=2)
context.map(process_cell, cellid_pattern)
# bad pixels
bp = np.zeros(offset.shape, np.uint32)
# offset related bad pixels
offset_mn = np.nanmedian(offset, axis=(0, 1))
offset_std = np.nanstd(offset, axis=(0, 1))
bp[(offset < offset_mn-thresholds_offset_sigma*offset_std) |
(offset > offset_mn+thresholds_offset_sigma*offset_std)] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value
bp[(offset < thresholds_offset_hard[0]) | (
offset > thresholds_offset_hard[1])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value
bp[~np.isfinite(offset)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
# noise related bad pixels
noise_mn = np.nanmedian(noise, axis=(0, 1))
noise_std = np.nanstd(noise, axis=(0, 1))
bp[(noise < noise_mn-thresholds_noise_sigma*noise_std) |
(noise > noise_mn+thresholds_noise_sigma*noise_std)] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
bp[(noise < thresholds_noise_hard[0]) | (
noise > thresholds_noise_hard[1])] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
bp[~np.isfinite(noise)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
idx = (cellid == cellid_pattern[0])
return offset, noise, channel, gg, bp, im[12, 12, idx], normal_test, cellid_pattern
```
%% Cell type:code id: tags:
``` python
offset_g = {}
noise_g = {}
badpix_g = {}
data_g = {}
ntest_g = {}
# Should be the same cell order for all modules & all gain stages
cellid_patterns_g = {}
inp = []
for gg, run_num in enumerate(run_nums):
run_path = Path(in_folder, f"r{run_num:04d}")
for i in modules:
inp.append((run_path, i, gg))
with multiprocessing.Pool(processes=parallel_num_procs) as pool:
results = pool.starmap(characterize_module, inp)
for ir, r in enumerate(results):
offset, noise, i, gg, bp, data, normal, cellid_pattern = r
if data is None:
warning(f"No data for module {i} of gain {gg}")
skip_plots = True
continue
qm = module_index_to_qm(i)
if qm not in offset_g:
offset_g[qm] = np.zeros(offset.shape[:3] + (3,))
print("Constant shape:", offset_g[qm].shape)
noise_g[qm] = np.zeros_like(offset_g[qm])
badpix_g[qm] = np.zeros_like(offset_g[qm], dtype=np.uint32)
data_g[qm] = np.full((ntrains, 3), np.nan)
ntest_g[qm] = np.zeros_like(offset_g[qm])
cellid_patterns_g[qm] = cellid_pattern
else:
if not np.array_equal(cellid_pattern, cellid_patterns_g[qm]):
raise ValueError("Inconsistent cell ID pattern between gain stages")
offset_g[qm][..., gg] = offset
noise_g[qm][..., gg] = noise
badpix_g[qm][..., gg] = bp
data_g[qm][:data.shape[0], gg] = data
ntest_g[qm][..., gg] = normal
hn, cn = np.histogram(data, bins=20)
print(f"{gain_names[gg]} gain, Module: {qm}. "
f"Number of processed trains per cell: {data.shape[0]}.")
```
%% Cell type:code id: tags:
``` python
if skip_plots:
import sys
print('Skipping plots')
sys.exit(0)
```
%% Cell type:code id: tags:
``` python
# Read report path and create file location tuple to add with the injection
proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
file_loc = raw_data_location_string(proposal, run_nums)
report = get_report(metadata_folder)
```
%% Cell type:code id: tags:
``` python
# TODO: add db_module when received from myMDC
# Create the modules dict of karabo_das and PDUs
qm_dict = {}
for i, k_da in zip(modules, karabo_da):
qm = module_index_to_qm(i)
qm_dict[qm] = {"karabo_da": k_da,
"db_module": ""}
```
%% Cell type:code id: tags:
``` python
from iCalibrationDB import OperatingCondition
def parallel_gain_parameter(condition):
if is_parallel_gain:
pgc = OperatingCondition()
pgc.name = 'Parallel gain'
pgc.description = ''
pgc.lower_deviation = 0
pgc.upper_deviation = 0
pgc.value = 1
condition.add_operating_condition(pgc)
```
%% Cell type:code id: tags:
``` python
# Retrieve existing constants for comparison
clist = ["Offset", "Noise", "BadPixelsDark"]
print('Retrieve pre-existing constants for comparison.')
old_const = {}
old_mdata = {}
for qm in offset_g.keys():
old_const[qm] = {}
old_mdata[qm] = {}
qm_db = qm_dict[qm]
karabo_da = qm_db["karabo_da"]
cellid_pattern = cellid_patterns_g[qm]
condition = Conditions.Dark.LPD(
memory_cells=max_cells,
bias_voltage=bias_voltage,
capacitor=capacitor_setting_s,
memory_cell_order=make_cell_order_condition(inject_cell_order, cellid_pattern),
)
parallel_gain_parameter(condition)
for const in clist:
constant = getattr(Constants.LPD, const)()
if not qm_db["db_module"]:
# This should be used in case of running notebook
# by a different method other than myMDC which already
# sends CalCat info.
qm_db["db_module"] = get_pdu_from_db(karabo_id, [karabo_da], constant,
condition, cal_db_interface,
snapshot_at=creation_time)[0]
data, mdata = get_from_db(karabo_id, karabo_da,
constant,
condition, None,
cal_db_interface,
creation_time=creation_time,
verbosity=2, timeout=cal_db_timeout)
old_const[qm][const] = data
if mdata is None or data is None:
old_mdata[qm][const] = {
"timestamp": "Not found",
"filepath": None,
"h5path": None
}
else:
timestamp = mdata.calibration_constant_version.begin_at.isoformat()
filepath = Path(
mdata.calibration_constant_version.hdf5path,
mdata.calibration_constant_version.filename
)
h5path = mdata.calibration_constant_version.h5path
old_mdata[qm][const] = {
"timestamp": timestamp,
"filepath": str(filepath),
"h5path": h5path
}
with open(f"{out_folder}/module_metadata_{qm}.yml","w") as fd:
yaml.safe_dump(
{
"module": qm,
"pdu": qm_db["db_module"],
"old-constants": old_mdata[qm]
}, fd)
```
%% Cell type:code id: tags:
``` python
res = {}
for i in modules:
qm = module_index_to_qm(i)
res[qm] = {'Offset': offset_g[qm],
'Noise': noise_g[qm],
'BadPixelsDark': badpix_g[qm]
}
```
%% Cell type:code id: tags:
``` python
# Save constants in the calibration DB
md = None
for qm in res:
karabo_da = qm_dict[qm]["karabo_da"]
db_module = qm_dict[qm]["db_module"]
mem_cell_order = make_cell_order_condition(
inject_cell_order, cellid_patterns_g[qm]
)
print("Memory cell order:", mem_cell_order)
# Do not store empty constants
# In case of 0 trains data_g is initiated with nans and never refilled.
if np.count_nonzero(~np.isnan(data_g[qm]))==0:
print(f"Constant ({qm}) would be empty, skipping saving")
continue
for const in res[qm]:
dconst = getattr(Constants.LPD, const)()
dconst.data = res[qm][const]
# set the operating condition
condition = Conditions.Dark.LPD(memory_cells=max_cells,
bias_voltage=bias_voltage,
capacitor=capacitor_setting_s,
memory_cell_order=mem_cell_order,
)
parallel_gain_parameter(condition)
if db_output:
md = send_to_db(db_module, karabo_id, dconst, condition,
file_loc, report_path=report,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout)
if local_output:
md = save_const_to_h5(db_module, karabo_id, dconst, condition,
dconst.data, file_loc, report, creation_time, out_folder)
print(f"Calibration constant {const} is stored locally.\n")
print("Constants parameter conditions are:\n")
print(f"• memory_cells: {max_cells}\n"
f"• bias_voltage: {bias_voltage}\n"
f"• capacitor: {capacitor_setting_s}\n"
f"• memory cell order: {mem_cell_order}\n"
f"• creation_time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\n")
```
%% Cell type:code id: tags:
``` python
show_processed_modules(
karabo_id,
constants=None,
mnames=[module_index_to_qm(i) for i in modules],
mode="position"
)
```
%% Cell type:markdown id: tags:
## Raw pedestal distribution ##
Distribution of a pedestal (ADUs) over trains for the pixel (12,12), memory cell 12. A median of the distribution is shown in yellow. A standard deviation is shown in red. The green line shows average over all pixels for a given memory cell and gain stage.
%% Cell type:code id: tags:
``` python
fig, grid = plt.subplots(3, 1, sharex="col", sharey="row", figsize=(10, 7))
fig.subplots_adjust(wspace=0, hspace=0)
for i in modules:
qm = module_index_to_qm(i)
if np.count_nonzero(~np.isnan(data_g[qm])) == 0:
break
for gain in range(3):
data = data_g[qm][:, gain]
offset = np.nanmedian(data)
noise = np.nanstd(data)
xrange = [np.nanmin(data_g[qm]), np.nanmax(data_g[qm])]
if xrange[1] == xrange[0]:
xrange = [0, xrange[0]+xrange[0]//2]
nbins = data_g[qm].shape[0]
else:
nbins = int(xrange[1] - xrange[0])
hn, cn = np.histogram(data, bins=nbins, range=xrange)
grid[gain].hist(data, range=xrange, bins=nbins)
grid[gain].plot([offset-noise, offset-noise], [0, np.nanmax(hn)],
linewidth=1.5, color='red',
label='1 $\sigma$ deviation')
grid[gain].plot([offset+noise, offset+noise],
[0, np.nanmax(hn)], linewidth=1.5, color='red')
grid[gain].plot([offset, offset], [0, 0],
linewidth=1.5, color='y', label='median')
grid[gain].plot([np.nanmedian(offset_g[qm][:, :, 12, gain]),
np.nanmedian(offset_g[qm][:, :, 12, gain])],
[0, np.nanmax(hn)], linewidth=1.5, color='green',
label='average over pixels')
grid[gain].set_xlim(xrange)
grid[gain].set_ylim(0, np.nanmax(hn)*1.1)
grid[gain].set_xlabel("Offset value [ADU]")
grid[gain].set_ylabel("# of occurance")
if gain == 0:
leg = grid[gain].legend(
loc='upper center', ncol=3,
bbox_to_anchor=(0.1, 0.25, 0.7, 1.0))
grid[gain].text(820, np.nanmax(hn)*0.4,
"{} gain".format(gain_names[gain]), fontsize=20)
a = plt.axes([.125, .1, 0.775, .8], frame_on=False)
a.patch.set_alpha(0.05)
a.set_xlim(xrange)
plt.plot([offset, offset], [0, 1], linewidth=1.5, color='y')
plt.xticks([])
plt.yticks([])
ypos = 0.9
x1pos = (np.nanmedian(data_g[qm][:, 0]) +
np.nanmedian(data_g[qm][:, 2]))/2.
x2pos = (np.nanmedian(data_g[qm][:, 2]) +
np.nanmedian(data_g[qm][:, 1]))/2.-10
plt.annotate("", xy=(np.nanmedian(data_g[qm][:, 0]), ypos), xycoords='data',
xytext=(np.nanmedian(data_g[qm][:, 2]), ypos), textcoords='data',
arrowprops=dict(arrowstyle="<->", connectionstyle="arc3"))
plt.annotate('{}'.format(np.nanmedian(data_g[qm][:, 0])-np.nanmedian(data_g[qm][:, 2])),
xy=(x1pos, ypos), xycoords='data', xytext=(5, 5), textcoords='offset points')
plt.annotate("", xy=(np.nanmedian(data_g[qm][:, 2]), ypos), xycoords='data',
xytext=(np.nanmedian(data_g[qm][:, 1]), ypos), textcoords='data',
arrowprops=dict(arrowstyle="<->", connectionstyle="arc3"))
plt.annotate('{}'.format(np.nanmedian(data_g[qm][:, 2])-np.nanmedian(data_g[qm][:, 1])),
xy=(x2pos, ypos), xycoords='data', xytext=(5, 5), textcoords='offset points')
plt.show()
```
%% Cell type:markdown id: tags:
## Normality test ##
Distributions of raw pedestal values have been tested if they are normally distributed. A normality test have been performed for each pixel and each memory cell. Plots below show histogram of p-Values and a 2D distribution for the memory cell 12.
%% Cell type:code id: tags:
``` python
# Loop over modules, constants
if not test_for_normality:
print('Normality test was not requested. Flag `test_for_normality` False')
else:
for i in modules:
qm = module_index_to_qm(i)
data = np.copy(ntest_g[qm][:,:,:,:])
data[badpix_g[qm][:,:,:,:]>0] = 1.01
hn,cn = np.histogram(data[:,:,:,0], bins=100)
d = [{'x': np.arange(100)*0.01+0.01,
'y': np.histogram(data[:,:,:,0], bins=100)[0],
'drawstyle': 'steps-pre',
'label' : 'High gain',
},
{'x': np.arange(100)*0.01+0.01,
'y': np.histogram(data[:,:,:,1], bins=100)[0],
'drawstyle': 'steps-pre',
'label' : 'Medium gain',
},
{'x': np.arange(100)*0.01+0.01,
'y': np.histogram(data[:,:,:,2], bins=100)[0],
'drawstyle': 'steps-pre',
'label' : 'Low gain',
},
]
fig = plt.figure(figsize=(15,15), tight_layout={'pad': 0.5, 'w_pad': 0.3})
for gain in range(3):
ax = fig.add_subplot(221+gain)
heatmapPlot(data[:,:,12,gain], add_panels=False, cmap='viridis', figsize=(10,10),
y_label='Rows', x_label='Columns',
lut_label='p-Value',
use_axis=ax,
title='p-Value for cell 12, {} gain'.format(gain_names[gain]) )
ax = fig.add_subplot(224)
_ = simplePlot(d, #aspect=1.6,
x_label = "p-Value".format(gain),
y_label="# of occurance",
use_axis=ax,
y_log=False, legend='outside-top-ncol3-frame', legend_pad=0.05, legend_size='5%')
ax.ticklabel_format(style='sci', axis='y', scilimits=(4,6))
```
%% Cell type:raw id: tags:
.. raw:: latex
\newpage
%% Cell type:markdown id: tags:
## Single-Cell Overviews ##
Single cell overviews allow to identify potential effects on all memory cells, e.g. on a sensor level. Additionally, they should serve as a first sanity check on expected behaviour, e.g. if structuring on the ASIC level is visible in the offsets, but otherwise no immediate artifacts are visible.
%% Cell type:code id: tags:
``` python
cell = 12
for gain in range(3):
display(
Markdown('### Cell-12 overview - {} gain'.format(gain_names[gain])))
fig = plt.figure(figsize=(18, 22) , tight_layout={'pad': 0.1, 'w_pad': 0.1})
for qm in res:
for iconst, const in enumerate(['Offset', 'Noise', 'BadPixelsDark']):
ax = fig.add_subplot(321+iconst)
data = res[qm][const][:, :, 12, gain]
vmax = 1.5 * np.nanmedian(res[qm][const][:, :, 12, gain])
title = const
label = '{} value [ADU]'.format(const)
title = '{} value'.format(const)
if const == 'BadPixelsDark':
vmax = 4
bpix_code = data.astype(np.float32)
bpix_code[bpix_code == 0] = np.nan
title = 'Bad pixel code'
label = title
cb_labels = ['1 {}'.format(BadPixels.NOISE_OUT_OF_THRESHOLD.name),
'2 {}'.format(BadPixels.OFFSET_NOISE_EVAL_ERROR.name),
'3 {}'.format(BadPixels.OFFSET_OUT_OF_THRESHOLD.name),
'4 {}'.format('MIXED')]
heatmapPlot(bpix_code, add_panels=False, cmap='viridis',
y_label='Rows', x_label='Columns',
lut_label='', vmax=vmax,
use_axis=ax, cb_ticklabels=cb_labels, cb_ticks = np.arange(4)+1,
title='{}'.format(title))
del bpix_code
else:
heatmapPlot(data, add_panels=False, cmap='viridis',
y_label='Rows', x_label='Columns',
lut_label=label, vmax=vmax,
use_axis=ax,
title='{}'.format(title))
for qm in res:
for iconst, const in enumerate(['Offset', 'Noise']):
data = res[qm][const]
dataBP = np.copy(data)
dataBP[res[qm]['BadPixelsDark'] > 0] = -1
x_ranges = [[0, 1500], [0, 40]]
hn, cn = np.histogram(
data[:, :, :, gain], bins=100, range=x_ranges[iconst])
hnBP, cnBP = np.histogram(dataBP[:, :, :, gain], bins=cn)
d = [{'x': cn[:-1],
'y': hn,
'drawstyle': 'steps-pre',
'label': 'All data',
},
{'x': cnBP[:-1],
'y': hnBP,
'drawstyle': 'steps-pre',
'label': 'Bad pixels masked',
},
]
ax = fig.add_subplot(325+iconst)
_ = simplePlot(d, figsize=(5, 7), aspect=1,
x_label="{} value [ADU]".format(const),
y_label="# of occurance",
title='', legend_pad=0.1, legend_size='10%',
use_axis=ax,
y_log=True, legend='outside-top-2col-frame')
plt.show()
```
%% Cell type:raw id: tags:
.. raw:: latex
\newpage
%% Cell type:markdown id: tags:
%% Cell type:code id: tags:
``` python
cols = {BadPixels.NOISE_OUT_OF_THRESHOLD.value: (BadPixels.NOISE_OUT_OF_THRESHOLD.name, '#FF000080'),
BadPixels.OFFSET_NOISE_EVAL_ERROR.value: (BadPixels.OFFSET_NOISE_EVAL_ERROR.name, '#0000FF80'),
BadPixels.OFFSET_OUT_OF_THRESHOLD.value: (BadPixels.OFFSET_OUT_OF_THRESHOLD.name, '#00FF0080'),
BadPixels.OFFSET_OUT_OF_THRESHOLD.value | BadPixels.NOISE_OUT_OF_THRESHOLD.value: ('MIXED', '#DD00DD80')}
if high_res_badpix_3d:
display(Markdown("""
## Global Bad Pixel Behaviour ##
The following plots shows the results of a bad pixel evaluation for all evaluated memory cells.
Cells are stacked in the Z-dimension, while pixels values in x/y are re-binned with a factor of 2.
This excludes single bad pixels present only in disconnected pixels.
Hence, any bad pixels spanning at least 4 pixels in the x/y-plane, or across at least two memory cells are indicated.
Colors encode the bad pixel type, or mixed type.
"""))
# Switch rebin to 1 for full resolution and
# no interpolation for badpixel values.
rebin = 2
for gain in range(3):
display(Markdown('### Bad pixel behaviour - {} gain ###'.format(gain_names[gain])))
for mod, data in badpix_g.items():
plot_badpix_3d(data[...,gain], cols, title='', rebin_fac=rebin)
ax = plt.gca()
leg = ax.get_legend()
leg.set(alpha=0.5)
plt.show()
```
%% Cell type:raw id: tags:
.. raw:: latex
\newpage
%% Cell type:markdown id: tags:
## Summary across tiles ##
Plots give an overview of calibration constants averaged across tiles. A bad pixel mask is applied. Constants are compared with pre-existing constants retrieved from the calibration database. Differences $\Delta$ between the old and new constants is shown.
%% Cell type:code id: tags:
``` python
time_summary = []
time_summary.append(f"The following pre-existing constants are used for comparison:")
for qm, qm_data in old_mdata.items():
time_summary.append(f"- Module {qm}")
for const, const_data in qm_data.items():
time_summary.append(f" - {const} created at {const_data['timestamp']}")
display(Markdown("\n".join(time_summary)))
```
%% Cell type:code id: tags:
``` python
# Loop over modules, constants
for qm in res:
for gain in range(3):
display(Markdown('### Summary across tiles - {} gain'.format(gain_names[gain])))
for const in res[qm]:
data = np.copy(res[qm][const][:, :, :, gain])
label = 'Fraction of bad pixels'
if const != 'BadPixelsDark':
data[badpix_g[qm][:, :, :, gain] > 0] = np.nan
label = '{} value [ADU]'.format(const)
else:
data[data>0] = 1.0
data = data.reshape(
int(data.shape[0] / 32),
32,
int(data.shape[1] / 128),
128,
data.shape[2])
data = np.nanmean(data, axis=(1, 3)).swapaxes(
0, 2).reshape(512, 16)
fig = plt.figure(figsize=(15, 6))
ax = fig.add_subplot(121)
_ = heatmapPlot(data[:510, :], add_panels=True,
y_label='Momery Cell ID', x_label='Tile ID',
lut_label=label, use_axis=ax,
panel_y_label=label, panel_x_label=label,
cmap='viridis', # cb_loc='right',cb_aspect=15,
x_ticklabels=np.arange(16)+1,
x_ticks=np.arange(16)+0.5)
if old_const[qm][const] is not None:
ax = fig.add_subplot(122)
dataold = np.copy(old_const[qm][const][:, :, :, gain])
label = '$\Delta$ {}'.format(label)
if const != 'BadPixelsDark':
if old_const[qm]['BadPixelsDark'] is not None:
dataold[old_const[qm]['BadPixelsDark'][:, :, :, gain] > 0] = np.nan
else:
dataold[:] = np.nan
else:
dataold[dataold>0]=1.0
dataold = dataold.reshape(
int(dataold.shape[0] / 32),
32,
int(dataold.shape[1] / 128),
128,
dataold.shape[2])
dataold = np.nanmean(dataold, axis=(
1, 3)).swapaxes(0, 2).reshape(512, 16)
dataold = dataold - data
_ = heatmapPlot(dataold[:510, :], add_panels=True,
y_label='Momery Cell ID', x_label='Tile ID',
lut_label=label, use_axis=ax,
panel_y_label=label, panel_x_label=label,
cmap='viridis', # cb_loc='right',cb_aspect=15,
x_ticklabels=np.arange(16)+1,
x_ticks=np.arange(16)+0.5)
plt.show()
```
%% Cell type:raw id: tags:
.. raw:: latex
\newpage
%% Cell type:markdown id: tags:
## Variation of offset and noise across Tiles and ASICs ##
The following plots show a standard deviation $\sigma$ of the calibration constant. The plot of standard deviations across tiles show pixels of one tile ($128 \times 32$). Value for each pixel shows a standard deviation across 16 tiles. The standard deviation across ASICs are shown overall tiles. The plot shows pixels of one ASIC ($16 \times 32$), where the value shows a standard deviation across all ACIS of the module.
%% Cell type:code id: tags:
``` python
# Loop over modules, constants
for qm in res:
for gain in range(3):
display(Markdown('### Variation of offset and noise across ASICs - {} gain'.format(gain_names[gain])))
fig = plt.figure(figsize=(15, 6))
for iconst, const in enumerate(['Offset', 'Noise']):
data = np.copy(res[qm][const][:, :, :, gain])
data[badpix_g[qm][:, :, :, gain] > 0] = np.nan
label = '$\sigma$ {} [ADU]'.format(const)
dataA = np.nanmean(data, axis=2) # average over cells
dataA = dataA.reshape(8, 32, 16, 16)
dataA = np.nanstd(dataA, axis=(0, 2)) # average across ASICs
ax = fig.add_subplot(121+iconst)
_ = heatmapPlot(dataA, add_panels=True,
y_label='rows', x_label='columns',
lut_label=label, use_axis=ax,
panel_y_label=label, panel_x_label=label,
cmap='viridis'
)
plt.show()
```
%% Cell type:code id: tags:
``` python
# Loop over modules, constants
for qm in res:
for gain in range(3):
display(Markdown('### Variation of offset and noise across tiles - {} gain'.format(gain_names[gain])))
fig = plt.figure(figsize=(15, 6))
for iconst, const in enumerate(['Offset', 'Noise']):
data = np.copy(res[qm][const][:, :, :, gain])
data[badpix_g[qm][:, :, :, gain] > 0] = np.nan
label = '$\sigma$ {} [ADU]'.format(const)
dataT = data.reshape(
int(data.shape[0] / 32),
32,
int(data.shape[1] / 128),
128,
data.shape[2])
dataT = np.nanstd(dataT, axis=(0, 2))
dataT = np.nanmean(dataT, axis=2)
ax = fig.add_subplot(121+iconst)
_ = heatmapPlot(dataT, add_panels=True,
y_label='rows', x_label='columns',
lut_label=label, use_axis=ax,
panel_y_label=label, panel_x_label=label,
cmap='viridis')
plt.show()
```
%% Cell type:raw id: tags:
.. raw:: latex
\newpage
%% Cell type:markdown id: tags:
## Aggregate values and per cell behaviour ##
The following tables and plots give an overview of statistical aggregates for each constant, as well as per-cell behavior, averaged across pixels.
%% Cell type:code id: tags:
``` python
# Loop over modules, constants
for qm in res:
for gain in range(3):
display(Markdown('### Mean over pixels - {} gain'.format(gain_names[gain])))
fig = plt.figure(figsize=(9,11))
for iconst, const in enumerate(res[qm]):
ax = fig.add_subplot(311+iconst)
data = res[qm][const][:,:,:510,gain]
if const == 'BadPixelsDark':
data[data>0] = 1.0
dataBP = np.copy(data)
dataBP[badpix_g[qm][:,:,:510,gain]>0] = -10
data = np.nanmean(data, axis=(0,1))
dataBP = np.nanmean(dataBP, axis=(0,1))
d = [{'y': data,
'x': np.arange(data.shape[0]),
'drawstyle': 'steps-mid',
'label' : 'All data'
}
]
if const != 'BadPixelsDark':
d.append({'y': dataBP,
'x': np.arange(data.shape[0]),
'drawstyle': 'steps-mid',
'label' : 'good pixels only'
})
y_title = "{} value [ADU]".format(const)
title = "{} value, {} gain".format(const, gain_names[gain])
else:
y_title = "Fraction of Bad Pixels"
title = "Fraction of Bad Pixels, {} gain".format(gain_names[gain])
data_min = np.min([data, dataBP])if const != 'BadPixelsDark' else np.min([data])
data_max = np.max([data[20:], dataBP[20:]])
data_dif = data_max - data_min
local_max = np.max([data[200:300], dataBP[200:300]])
frac = 0.35
new_max = (local_max - data_min*(1-frac))/frac
new_max = np.max([data_max, new_max])
_ = simplePlot(d, figsize=(10,10), aspect=2, xrange=(-12, 510),
x_label = 'Memory Cell ID',
y_label=y_title, use_axis=ax,
title=title,
title_position=[0.5, 1.15],
inset='xy-coord-right', inset_x_range=(0,20), inset_indicated=True,
inset_labeled=True, inset_coord=[0.2,0.5,0.6,0.95],
inset_lw = 1.0, y_range = [data_min-data_dif*0.05, new_max+data_dif*0.05],
y_log=False, legend='outside-top-ncol2-frame', legend_size='18%',
legend_pad=0.00)
plt.tight_layout(pad=1.08, h_pad=0.35)
plt.show()
```
%% Cell type:raw id: tags:
.. raw:: latex
\newpage
%% Cell type:markdown id: tags:
## Summary tables ##
The following tables show summary information for the evaluated module. Values for currently evaluated constants are compared with values for pre-existing constants retrieved from the calibration database.
%% Cell type:code id: tags:
``` python
table = []
bits = [
BadPixels.NOISE_OUT_OF_THRESHOLD,
BadPixels.OFFSET_OUT_OF_THRESHOLD,
BadPixels.OFFSET_NOISE_EVAL_ERROR
]
for qm in res:
for gain in range(3):
l_data = []
l_data_old = []
data = np.copy(res[qm]['BadPixelsDark'][:,:,:,gain])
l_data.append(len(data[data>0].flatten()))
for bit in bits:
l_data.append(np.count_nonzero(badpix_g[qm][:,:,:,gain] & bit.value))
if old_const[qm]['BadPixelsDark'] is not None:
old_const[qm]['BadPixelsDark'] = old_const[qm]['BadPixelsDark'].astype(np.uint32)
dataold = np.copy(old_const[qm]['BadPixelsDark'][:, :, :, gain])
l_data_old.append(len(dataold[dataold>0].flatten()))
for bit in bits:
l_data_old.append(np.count_nonzero(old_const[qm]['BadPixelsDark'][:, :, :, gain] & bit.value))
l_data_name = ['All bad pixels', 'NOISE_OUT_OF_THRESHOLD',
'OFFSET_OUT_OF_THRESHOLD', 'OFFSET_NOISE_EVAL_ERROR']
l_threshold = ['', f'{thresholds_noise_sigma}', f'{thresholds_offset_sigma}',
f'{thresholds_offset_hard}/{thresholds_noise_hard}']
for i in range(len(l_data)):
line = [f'{l_data_name[i]}, gain {gain_names[gain]}', l_threshold[i], l_data[i]]
if old_const[qm]['BadPixelsDark'] is not None:
line += [l_data_old[i]]
else:
line += ['-']
table.append(line)
table.append(['', '', '', ''])
display(Markdown('''
### Number of bad pixels ###
One pixel can be bad for different reasons, therefore, the sum of all types of bad pixels can be more than the number of all bad pixels.
'''))
if len(table)>0:
md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
headers=["Pixel type", "Threshold",
"New constant", "Old constant"])))
```
%% Cell type:code id: tags:
``` python
header = ['Parameter',
"New constant", "Old constant ",
"New constant", "Old constant ",
"New constant", "Old constant "]
for const in ['Offset', 'Noise']:
table = [['','High gain', 'High gain', 'Medium gain', 'Medium gain', 'Low gain', 'Low gain']]
for qm in res:
data = np.copy(res[qm][const])
data[res[qm]['BadPixelsDark']>0] = np.nan
if old_const[qm][const] is not None and old_const[qm]['BadPixelsDark'] is not None :
dataold = np.copy(old_const[qm][const])
dataold[old_const[qm]['BadPixelsDark']>0] = np.nan
f_list = [np.nanmedian, np.nanmean, np.nanstd, np.nanmin, np.nanmax]
n_list = ['Median', 'Mean', 'Std', 'Min', 'Max']
for i, f in enumerate(f_list):
line = [n_list[i]]
for gain in range(3):
line.append('{:6.1f}'.format(f(data[...,gain])))
if old_const[qm][const] is not None and old_const[qm]['BadPixelsDark'] is not None:
line.append('{:6.1f}'.format(f(dataold[...,gain])))
else:
line.append('-')
table.append(line)
display(Markdown('### {} [ADU], good pixels only ###'.format(const)))
md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=header)))
```

Legal Notice | Contacts