Skip to content
Snippets Groups Projects
Commit f80914a6 authored by Mikhail Karnevskiy's avatar Mikhail Karnevskiy
Browse files

Fix array size

parent 9d5b8ab0
No related branches found
No related tags found
1 merge request!111Fixes inspired by deployment
%% Cell type:markdown id: tags:
# 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.
**The offset** ($O$) is defined as the median ($M$) of the dark signal ($Ds$) over trains ($t$) for a given pixel ($x,y$) and memory cell ($c$).
**The noise** $N$ is the standard deviation $\sigma$ of the dark signal.
$$ O_{x,y,c} = M(Ds)_{t} ,\,\,\,\,\,\, N_{x,y,c} = \sigma(Ds)_{t}$$
**The bad pixel** mask is encoded as a bit mask.
**"OFFSET_OUT_OF_THRESHOLD":**
Offset outside of bounds:
$$M(O)_{x,y} - \sigma(O)_{x,y} * \mathrm{thresholds\_offset\_sigma} < O < M(O)_{x,y} + \sigma(O)_{x,y} * \mathrm{thresholds\_offset\_sigma} $$
or offset outside of hard limits
$$ \mathrm{thresholds\_offset\_hard}_\mathrm{low} < O < \mathrm{thresholds\_offset\_hard}_\mathrm{high} $$
**"NOISE_OUT_OF_THRESHOLD":**
Noise outside of bounds:
$$M(N)_{x,y} - \sigma(N)_{x,y} * \mathrm{thresholds\_noise\_sigma} < N < M(N)_{x,y} + \sigma(N)_{x,y} * \mathrm{thresholds\_noise\_sigma} $$
or noise outside of hard limits
$$\mathrm{thresholds\_noise\_hard}_\mathrm{low} < N < \mathrm{thresholds\_noise\_hard}_\mathrm{high} $$
**"OFFSET_NOISE_EVAL_ERROR":**
Offset and Noise both not $nan$ values
Values: $\mathrm{thresholds\_offset\_sigma}$, $\mathrm{thresholds\_offset\_hard}$, $\mathrm{thresholds\_noise\_sigma}$, $\mathrm{thresholds\_noise\_hard}$ are given as parameters.
%% Cell type:code id: tags:
``` python
cluster_profile = "noDB" # The ipcluster profile to use
in_folder = "/gpfs/exfel/exp/FXE/201930/p900063/raw" # path to input data, required
in_folder = "/gpfs/exfel/exp/FXE/201931/p900088/raw" # path to input data, required
out_folder = "/gpfs/exfel/data/scratch/karnem/LPD/" # path to output to, required
sequences = [0] # sequence files to evaluate
modules = [-1] # list of modules to evaluate, RANGE ALLOWED
capacitor_setting = 5 # capacitor_setting for which data was taken, required
run_high = 358 # run number in which high gain data was recorded, required
run_med = 359 # run number in which medium gain data was recorded, required
run_low = 360 # run number in which low gain data was recorded, required
run_high = 112 # run number in which high gain data was recorded, required
run_med = 113 # run number in which medium gain data was recorded, required
run_low = 114 # run number in which low gain data was recorded, required
mem_cells = 512 # number of memory cells used
local_output = True # output constants locally
db_output = True # output constants to database
bias_voltage = 250 # detector bias voltage
cal_db_interface = "tcp://max-exfl016:8015#8025" # the database interface to use
cal_db_timeout = 300000 # timeout on caldb requests"
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
skip_first_ntrains = 10 # Number of first trains to skip
not_use_dir_creation_date = False # do not use the creation date of the directory for database time derivation
instrument = "FXE" # instrument name
ntrains = 300 # number of trains to use
high_res_badpix_3d = False # plot bad-pixel summary in high resolution
```
%% Cell type:code id: tags:
``` python
from collections import OrderedDict
import copy
from datetime import datetime
from functools import partial
import os
import warnings
warnings.filterwarnings('ignore')
import dateutil.parser
import h5py
from ipyparallel import Client
from IPython.display import display, Markdown, Latex
import numpy as np
import matplotlib
matplotlib.use("agg")
import matplotlib.patches as patches
import matplotlib.pyplot as plt
%matplotlib inline
import tabulate
from XFELDetAna.plotting.heatmap import heatmapPlot
from XFELDetAna.plotting.simpleplot import simplePlot
from iCalibrationDB import (ConstantMetaData, Constants,
Conditions, Detectors,
Versions)
from cal_tools.tools import (gain_map_files, parse_runs,
run_prop_seq_from_path,
get_notebook_name,
get_dir_creation_date, get_from_db)
get_dir_creation_date, get_from_db,
get_random_db_interface)
from cal_tools.influx import InfluxLogger
from cal_tools.enums import BadPixels
from cal_tools.plotting import (show_overview, plot_badpix_3d,
create_constant_overview)
```
%% Cell type:code id: tags:
``` python
client = Client(profile=cluster_profile)
view = client[:]
view.use_dill()
gains = np.arange(3)
max_cells = mem_cells
cells = np.arange(max_cells)
gain_names = ['High', 'Medium', 'Low']
if modules[0] == -1:
modules = list(range(16))
gain_runs = OrderedDict()
if capacitor_setting == 5:
gain_runs["high_5pf"] = "r{:04d}".format(run_high)
gain_runs["med_5pf"] = "r{:04d}".format(run_med)
gain_runs["low_5pf"] = "r{:04d}".format(run_low)
elif capacitor_setting == 50:
gain_runs["high_50pf"] = "r{:04d}".format(run_high)
gain_runs["med_50pf"] = "r{:04d}".format(run_med)
gain_runs["low_50pf"] = "r{:04d}".format(run_low)
capacitor_settings = [capacitor_setting]
capacitor_settings = ['{}pf'.format(c) for c in capacitor_settings]
QUADRANTS = 4
MODULES_PER_QUAD = 4
DET_FILE_INSET = "LPD"
creation_time = None
if not not_use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run_high)
print("Using {} as creation time".format(creation_time))
run, prop, seq = run_prop_seq_from_path(in_folder)
logger = InfluxLogger(detector="LPD", instrument=instrument, mem_cells=mem_cells,
notebook=get_notebook_name(), proposal=prop)
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("Runs: {}, {}, {}".format(run_high, run_med, run_low))
print("Sequences: {}".format(sequences))
print("Using DB: {}".format(db_output))
print("Input: {}".format(in_folder))
print("Output: {}".format(out_folder))
print("Bias voltage: {}V".format(bias_voltage))
```
%% Cell type:code id: tags:
``` python
# set everything up filewise
if not os.path.exists(out_folder):
os.makedirs(out_folder)
gmf = gain_map_files(in_folder, gain_runs, sequences, DET_FILE_INSET, QUADRANTS, MODULES_PER_QUAD)
gain_mapped_files, total_sequences, total_file_size = gmf
```
%% Cell type:markdown id: tags:
## Data processing
%% Cell type:code id: tags:
``` python
# the actual characterization - to not eded this without consultation
def characterize_module(cells, bp_thresh, skip_first_ntrains, ntrains, inp):
import numpy as np
import copy
import h5py
from cal_tools.enums import BadPixels
import scipy.stats
def splitOffGainLPD(d):
msk = np.zeros(d.shape, np.uint16)
msk[...] = 0b0000111111111111
data = np.bitwise_and(d, msk)
msk[...] = 0b0011000000000000
gain = np.bitwise_and(d, msk)//4096
gain[gain > 2] = 2
return data, gain
filename, filename_out, channel = inp
thresholds_offset_hard, thresholds_offset_sigma, thresholds_noise_hard, thresholds_noise_sigma = bp_thresh
infile = h5py.File(filename, "r", driver="core")
im = np.array(infile["/INSTRUMENT/FXE_DET_LPD1M-1/DET/{}CH0:xtdf/image/data".format(
channel)][skip_first_ntrains*cells:skip_first_ntrains*cells+ntrains*cells, ...])
cellid = np.squeeze(np.array(infile["/INSTRUMENT/FXE_DET_LPD1M-1/DET/{}CH0:xtdf/image/cellId".format(
channel)][skip_first_ntrains*cells:skip_first_ntrains*cells+ntrains*cells, ...]))
infile.close()
im, g = splitOffGainLPD(im[:, 0, ...])
im = im.astype(np.float32)
im = np.rollaxis(im, 2)
im = np.rollaxis(im, 2, 1)
offset = np.zeros((im.shape[0], im.shape[1], cells))
noise = np.zeros((im.shape[0], im.shape[1], cells))
normal_test = np.zeros((im.shape[0], im.shape[1], cells))
for cc in range(cells):
idx = cellid == cc
if np.any(idx):
offset[..., cc] = np.median(im[:, :, idx], axis=2)
noise[..., cc] = np.std(im[:, :, idx], axis=2)
_, normal_test[..., cc] = scipy.stats.normaltest(
im[:, :, idx], axis=2)
# 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 == 12
return offset, noise, channel, bp, im[12, 12, idx], normal_test
offset_g = OrderedDict()
noise_g = OrderedDict()
badpix_g = OrderedDict()
data_g = OrderedDict()
ntest_g = OrderedDict()
gg = 0
old_cap = None
start = datetime.now()
for gain, mapped_files in gain_mapped_files.items():
cap = gain.split("_")[1]
if cap != old_cap:
gg = 0
old_cap = cap
offset_g[cap] = OrderedDict()
noise_g[cap] = OrderedDict()
badpix_g[cap] = OrderedDict()
data_g[cap] = OrderedDict()
ntest_g[cap] = OrderedDict()
dones = []
inp = []
for i in modules:
qm = "Q{}M{}".format(i//4 + 1, i % 4 + 1)
if qm in mapped_files and not mapped_files[qm].empty():
fname_in = mapped_files[qm].get()
dones.append(mapped_files[qm].empty())
else:
continue
fout = os.path.abspath(
"{}/{}".format(out_folder, (os.path.split(fname_in)[-1]).replace("RAW", "CORR")))
print("Process file: ", fout)
inp.append((fname_in, fout, i))
first = False
p = partial(characterize_module, max_cells,
(thresholds_offset_hard, thresholds_offset_sigma,
thresholds_noise_hard, thresholds_noise_sigma),
skip_first_ntrains, ntrains)
results = view.map_sync(p, inp)
for r in results:
offset, noise, i, bp, data, normal = r
qm = "Q{}M{}".format(i//4 + 1, i % 4 + 1)
if qm not in offset_g[cap]:
offset_g[cap][qm] = np.zeros(
(offset.shape[0], offset.shape[1], offset.shape[2], 3))
noise_g[cap][qm] = np.zeros_like(offset_g[cap][qm])
badpix_g[cap][qm] = np.zeros_like(offset_g[cap][qm])
data_g[cap][qm] = np.zeros((data.shape[0], 3))
data_g[cap][qm] = np.zeros((ntrains, 3))
ntest_g[cap][qm] = np.zeros_like(offset_g[cap][qm])
offset_g[cap][qm][..., gg] = offset
noise_g[cap][qm][..., gg] = noise
badpix_g[cap][qm][..., gg] = bp
data_g[cap][qm][..., gg] = data
data_g[cap][qm][:data.shape[0], gg] = data
ntest_g[cap][qm][..., gg] = normal
hn, cn = np.histogram(data, bins=20)
print("{} gain. Module: {}. Number of processed trains per cell: {}.\n".format(
gain_names[gg], qm, data.shape[0]))
gg += 1
plt.show()
duration = (datetime.now()-start).total_seconds()
logger.runtime_summary_entry(success=True, runtime=duration,
total_sequences=total_sequences,
filesize=total_file_size)
logger.send()
```
%% Cell type:code id: tags:
``` python
# save everything to file.
if local_output:
for cap in capacitor_settings:
runs = [v for k, v in gain_runs.items() if cap in k]
ofile = "{}/lpd_offset_store_{}_{}_{}.h5".format(out_folder,
"_".join(runs),
cap,
"_".join([str(m) for m in modules]))
store_file = h5py.File(ofile, "w")
for qm in offset_g[cap].keys():
store_file["{}/Offset/0/data".format(qm)] = offset_g[cap][qm]
store_file["{}/Noise/0/data".format(qm)] = noise_g[cap][qm]
store_file["{}/BadPixelsDark/0/data".format(qm)
] = badpix_g[cap][qm]
store_file.close()
print('Constants are stored to {}'.format(ofile))
```
%% Cell type:code id: tags:
``` python
# Retrieve existing constants for comparison
clist = ["Offset", "Noise", "BadPixelsDark"]
old_const = {}
old_mdata = {}
print('Retrieve pre-existing constants for comparison.')
for cap in capacitor_settings:
for qm in offset_g[cap].keys():
for const in clist:
condition = Conditions.Dark.LPD(memory_cells=max_cells, bias_voltage=bias_voltage,
capacitor=cap)
data, mdata = get_from_db(getattr(Detectors.LPD1M1, qm),
getattr(Constants.LPD, const)(),
condition,
None,
cal_db_interface, creation_time=creation_time,
verbosity=2, timeout=cal_db_timeout)
old_const[const] = data
if mdata is not None and data is not None:
time = mdata.calibration_constant_version.begin_at
old_mdata[const] = time.isoformat()
else:
old_mdata[const] = "Not found"
```
%% Cell type:code id: tags:
``` python
res = OrderedDict()
for cap in capacitor_settings:
res[cap] = OrderedDict()
for i in modules:
qm = "Q{}M{}".format(i//4+1, i % 4+1)
res[cap][qm] = {'Offset': offset_g[cap][qm],
'Noise': noise_g[cap][qm],
'BadPixelsDark': badpix_g[cap][qm]
}
```
%% Cell type:code id: tags:
``` python
# Save constants in the calibration DB
if db_output:
for cap in capacitor_settings:
for qm in res[cap]:
for const in res[cap][qm]:
metadata = ConstantMetaData()
dconst = getattr(Constants.LPD, const)()
dconst.data = res[cap][qm][const]
metadata.calibration_constant = dconst
# set the operating condition
condition = Conditions.Dark.LPD(memory_cells=max_cells,
bias_voltage=bias_voltage,
capacitor=cap)
device = getattr(Detectors.LPD1M1, qm)
if device:
metadata.detector_condition = condition
# specify the a version for this constant
if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device)
else:
metadata.calibration_constant_version = Versions.Timespan(device=device,
start=creation_time)
metadata.send(cal_db_interface, timeout=cal_db_timeout)
msg = 'Const {} for module {} was injected to the calibration DB. Begin at: {}'
print(msg.format(const, qm,
metadata.calibration_constant_version.begin_at))
```
%% Cell type:code id: tags:
``` python
qm = "Q{}M{}".format(modules[0]//4+1, modules[0] % 4+1)
display(Markdown('## Position of the module {}, it tiles and ASICs of tile ##'.format(qm)))
fig, ax = plt.subplots(1, figsize=(10, 10))
ax.set_axis_off()
ax.set_xlim(0, 97)
ax.set_ylim(0, 97)
q_poses = np.array([[51, 47], [47, 1], [1, 5], [5, 51]])
m_poses = np.array([[22.5, 20.5], [22.5, 0.5], [0.5, 0.5], [0.5, 20.5]])
for iq, q_pos in enumerate(q_poses):
ax.add_patch(patches.Rectangle(q_pos, 45, 45, linewidth=2, edgecolor='r',
facecolor='y', fill=True))
ax.text(q_pos[0]+20, q_pos[1]+41.5, 'Q{}'.format(iq+1), fontsize=22)
for im, m_pos in enumerate(m_poses):
ax.add_patch(patches.Rectangle(q_pos+m_pos, 22, 20, linewidth=3, edgecolor='r',
facecolor='g', fill=True))
if iq*4+im == modules[0]:
for a_posx in range(2):
for a_posy in range(8):
a_pos = np.array([a_posx*11., a_posy*20/8.])
pos = q_pos+m_pos+a_pos
ax.add_patch(patches.Rectangle(q_pos+m_pos+a_pos, 11, 20/8.,
linewidth=1, edgecolor='black',
facecolor='r', fill=True))
if a_posx == 0:
label = str(a_posy+9)
else:
label = str(-a_posy+(a_posx*8))
ax.text(pos[0]+4, pos[1]+0.3, label, fontsize=14)
else:
pos = q_pos+m_pos+np.array([5, 8])
ax.text(pos[0], pos[1], 'Q{}M{}'.format(
iq+1, im+1), fontsize=22, color='y')
ax.add_patch(patches.Rectangle([65, 93], 30, 4, linewidth=1, edgecolor='black',
facecolor='r', fill=True))
ax.text(52, 94, 'ASICs:', fontsize=22, color='black')
for i_pos in range(8):
pos = np.array([65, 93]) + np.array([i_pos*30/8.+0.3, 0.3])
ax.add_patch(patches.Rectangle(pos, 24/8., 3.4, linewidth=1, edgecolor='black',
facecolor='deepskyblue', fill=True))
ax.text(pos[0]+0.5, pos[1]+0.5, '{}'.format(i_pos + 1),
fontsize=18, color='black')
```
%% 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 cap in capacitor_settings:
for i in modules:
qm = "Q{}M{}".format(i//4+1, i % 4+1)
if data_g[cap][qm].shape[0] == 0:
break
for gain in range(3):
data = data_g[cap][qm][:, gain]
offset = np.nanmedian(data)
noise = np.nanstd(data)
xrange = [np.nanmin(data_g[cap][qm]), np.nanmax(data_g[cap][qm])]
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[cap][qm][:, :, 12, gain]),
np.nanmedian(offset_g[cap][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='outside-top', 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[cap][qm][:, 0]) +
np.nanmedian(data_g[cap][qm][:, 2]))/2.
x2pos = (np.nanmedian(data_g[cap][qm][:, 2]) +
np.nanmedian(data_g[cap][qm][:, 1]))/2.-10
plt.annotate("", xy=(np.nanmedian(data_g[cap][qm][:, 0]), ypos), xycoords='data',
xytext=(np.nanmedian(data_g[cap][qm][:, 2]), ypos), textcoords='data',
arrowprops=dict(arrowstyle="<->", connectionstyle="arc3"))
plt.annotate('{}'.format(np.nanmedian(data_g[cap][qm][:, 0])-np.nanmedian(data_g[cap][qm][:, 2])),
xy=(x1pos, ypos), xycoords='data', xytext=(5, 5), textcoords='offset points')
plt.annotate("", xy=(np.nanmedian(data_g[cap][qm][:, 2]), ypos), xycoords='data',
xytext=(np.nanmedian(data_g[cap][qm][:, 1]), ypos), textcoords='data',
arrowprops=dict(arrowstyle="<->", connectionstyle="arc3"))
plt.annotate('{}'.format(np.nanmedian(data_g[cap][qm][:, 2])-np.nanmedian(data_g[cap][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 capacitor settings, modules, constants
for cap in capacitor_settings:
for i in modules:
qm = "Q{}M{}".format(i//4+1, i%4+1)
data = np.copy(ntest_g[cap][qm][:,:,:,:])
data[badpix_g[cap][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 cap in capacitor_settings:
for gain in range(3):
display(
Markdown('### Cell-12 overview - {} gain'.format(gain_names[gain])))
fig = plt.figure(figsize=(15, 12) , tight_layout={'pad': 0.1, 'w_pad': 0.1})
for qm in res[cap]:
for iconst, const in enumerate(['Offset', 'Noise', 'BadPixelsDark']):
ax = fig.add_subplot(221+iconst)
data = res[cap][qm][const][:, :, 12, gain]
vmax = 1.5 * np.nanmedian(res[cap][qm][const][:, :, 12, gain])
title = const
label = '{} value [ADU]'.format(const)
title = '{} value'.format(const)
if const == 'BadPixelsDark':
vmax = 4
data[data == 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(data, 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))
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))
#show_overview(res[cap], cell, gain, out_folder=out_folder, infix="_".join(gain_runs.values()))
fig = plt.figure(figsize=(10, 5))
for qm in res[cap]:
for iconst, const in enumerate(['Offset', 'Noise']):
data = res[cap][qm][const]
dataBP = np.copy(data)
dataBP[res[cap][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(121+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:
## 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.
%% 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')}
rebin = 8 if not high_res_badpix_3d else 2
for gain in range(3):
display(Markdown('### Bad pixel behaviour - {} gain ###'.format(gain_names[gain])))
for cap in capacitor_settings:
for mod, data in badpix_g[cap].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
display(Markdown('The following pre-existing constants are used for comparison: \n'))
for key in old_mdata:
display(Markdown('**{}** at {}'.format(key, old_mdata[key])))
```
%% Cell type:code id: tags:
``` python
# Loop over capacitor settings, modules, constants
for cap in res:
for qm in res[cap]:
for gain in range(3):
display(Markdown('### Summary across tiles - {} gain'.format(gain_names[gain])))
for const in res[cap][qm]:
data = np.copy(res[cap][qm][const][:, :, :, gain])
label = 'Fraction of bad pixels'
if const != 'BadPixelsDark':
data[badpix_g[cap][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[const] is not None:
ax = fig.add_subplot(122)
dataold = np.copy(old_const[const][:, :, :, gain])
label = '$\Delta$ {}'.format(label)
if const != 'BadPixelsDark':
if old_const['BadPixelsDark'] is not None:
dataold[old_const['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 capacitor settings, modules, constants
for cap in res:
for qm in res[cap]:
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[cap][qm][const][:, :, :, gain])
data[badpix_g[cap][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 capacitor settings, modules, constants
for cap in res:
for qm in res[cap]:
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[cap][qm][const][:, :, :, gain])
data[badpix_g[cap][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 capacitor settings, modules, constants
for cap in res:
for qm in res[cap]:
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[cap][qm]):
ax = fig.add_subplot(311+iconst)
data = res[cap][qm][const][:,:,:510,gain]
if const == 'BadPixelsDark':
data[data>0] = 1.0
dataBP = np.copy(data)
dataBP[badpix_g[cap][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 = []
for cap in res:
for qm in res[cap]:
for gain in range(3):
l_data = []
l_data_old = []
data = np.copy(res[cap][qm]['BadPixelsDark'][:,:,:,gain])
datau32 = data.astype(np.uint32)
l_data.append(data)
l_data.append(datau32 - np.bitwise_or(datau32,BadPixels.NOISE_OUT_OF_THRESHOLD.value))
l_data.append(datau32 - np.bitwise_or(datau32,BadPixels.OFFSET_OUT_OF_THRESHOLD.value))
l_data.append(datau32 - np.bitwise_or(datau32,BadPixels.OFFSET_NOISE_EVAL_ERROR.value))
if old_const['BadPixelsDark'] is not None:
dataold = np.copy(old_const['BadPixelsDark'][:, :, :, gain])
datau32old = dataold.astype(np.uint32)
l_data_old.append(dataold)
l_data_old.append(datau32old - np.bitwise_or(datau32old,BadPixels.NOISE_OUT_OF_THRESHOLD.value))
l_data_old.append(datau32old - np.bitwise_or(datau32old,BadPixels.OFFSET_OUT_OF_THRESHOLD.value))
l_data_old.append(datau32old - np.bitwise_or(datau32old,BadPixels.OFFSET_NOISE_EVAL_ERROR.value))
l_data_name = ['All bad pixels', 'NOISE_OUT_OF_THRESHOLD',
'OFFSET_OUT_OF_THRESHOLD', 'OFFSET_NOISE_EVAL_ERROR']
l_threshold = ['', '{}'.format(thresholds_noise_sigma), '{}'.format(thresholds_offset_sigma),
'{}/{}'.format(thresholds_offset_hard, thresholds_noise_hard)]
for i in range(len(l_data)):
line = ['{}, gain {}'.format(l_data_name[i], gain_names[gain]),
l_threshold[i],
len(l_data[i][l_data[i]>0].flatten())
]
if old_const['BadPixelsDark'] is not None:
line += [len(l_data_old[i][l_data_old[i]>0].flatten())]
else:
line += ['-']
table.append(line)
display(Markdown('### Number of bad pixels ###'.format(qm)))
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 cap in res:
for qm in res[cap]:
data = np.copy(res[cap][qm][const])
data[res[cap][qm]['BadPixelsDark']>0] = np.nan
if old_const[const] is not None and old_const['BadPixelsDark'] is not None :
dataold = np.copy(old_const[const])
dataold[old_const['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[const] is not None and old_const['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)))
```
......
......@@ -10,9 +10,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"cluster_profile = \"noDB\" # The ipcluster profile to use\n",
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment