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

Introduce fixes with respect to comments

parent 89080e10
No related branches found
No related tags found
1 merge request!92Feat: update plotting notebooks
Hide whitespace changes
Inline Side-by-side
cal_tools/cal_tools/ana_tools.py
+ 25
− 25
View file @ bff91e27
...@@ -2,6 +2,7 @@ import h5py ...@@ -2,6 +2,7 @@ import h5py
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import datetime import datetime
import dateutil.parser
import glob import glob
...@@ -43,16 +44,16 @@ def load_data_from_hdf5(filelist): ...@@ -43,16 +44,16 @@ def load_data_from_hdf5(filelist):
if mKey not in data[cKey]: if mKey not in data[cKey]:
data[cKey][mKey] = {} data[cKey][mKey] = {}
# Loop over module data # Loop over module data
for tKey in f.get(cKey + '/' + mKey): for tKey in f.get('{}/{}'.format(cKey, mKey)):
if tKey == 'ctime': if tKey == 'ctime':
if "ctime" not in data[cKey][mKey]: ctime = data[cKey][mKey].get("ctime", [])
data[cKey][mKey]["ctime"] = []
for timeKey in f.get("/".join((cKey, mKey, tKey))): for timeKey in f.get("/".join((cKey, mKey, tKey))):
path = "/".join((cKey, mKey, tKey, timeKey)) path = "/".join((cKey, mKey, tKey, timeKey))
value = f.get(path).value value = f.get(path).value
value = datetime.datetime.fromtimestamp(value) value = dateutil.parser.parse(value)
data[cKey][mKey]["ctime"].append(value) ctime.append(value)
data[cKey][mKey]["ctime"] = ctime
continue continue
# Load ndarray if they are there # Load ndarray if they are there
...@@ -63,7 +64,7 @@ def load_data_from_hdf5(filelist): ...@@ -63,7 +64,7 @@ def load_data_from_hdf5(filelist):
continue continue
# Loop over stored data # Loop over stored data
for dKay in f.get( "/".join((cKey, mKey, tKey))): for dKay in f.get("/".join((cKey, mKey, tKey))):
# Loop over metadata # Loop over metadata
if dKay == "mdata": if dKay == "mdata":
...@@ -74,19 +75,19 @@ def load_data_from_hdf5(filelist): ...@@ -74,19 +75,19 @@ def load_data_from_hdf5(filelist):
"/".join( "/".join(
(cKey, mKey, tKey, 'mdata', (cKey, mKey, tKey, 'mdata',
mdKey))).value mdKey))).value
if "mdata" not in data[cKey][mKey]: mdata_l = data[cKey][mKey].get("mdata", [])
data[cKey][mKey]["mdata"] = [] mdata_l.append(mdata)
data[cKey][mKey]["mdata"].append(mdata) data[cKey][mKey]["mdata"] = mdata_l
continue continue
if dKay not in data[cKey][mKey]: if dKay not in data[cKey][mKey]:
data[cKey][mKey][dKay] = [] data[cKey][mKey][dKay] = []
value = f.get( value = f.get(
"/".join((cKey, mKey, tKey, dKay))).value "/".join((cKey, mKey, tKey, dKay))).value
if dKay == "ctime": if dKay == "ctime":
value = datetime.datetime.fromtimestamp(value) value = dateutil.parser.parse(value)
data[cKey][mKey][dKay].append(value) data[cKey][mKey][dKay].append(value)
return data return data
...@@ -106,7 +107,7 @@ def recursively_save_dict_contents_to_group(h5file, path, dic): ...@@ -106,7 +107,7 @@ def recursively_save_dict_contents_to_group(h5file, path, dic):
# Prepare datatime to save # Prepare datatime to save
if isinstance(item, datetime.datetime): if isinstance(item, datetime.datetime):
item = item.timestamp() item = item.isoformat()
# Iterate over dictionary # Iterate over dictionary
if isinstance(item, dict): if isinstance(item, dict):
...@@ -119,7 +120,7 @@ def recursively_save_dict_contents_to_group(h5file, path, dic): ...@@ -119,7 +120,7 @@ def recursively_save_dict_contents_to_group(h5file, path, dic):
if isinstance(x, dict): if isinstance(x, dict):
recursively_save_dict_contents_to_group(h5file, newpath, x) recursively_save_dict_contents_to_group(h5file, newpath, x)
elif isinstance(x, datetime.datetime): elif isinstance(x, datetime.datetime):
h5file[newpath] = x.timestamp() h5file[newpath] = x.isoformat()
else: else:
raise ValueError('Cannot save a list of %s ' % type(item)) raise ValueError('Cannot save a list of %s ' % type(item))
...@@ -302,11 +303,11 @@ def multi_intersect(a, b): ...@@ -302,11 +303,11 @@ def multi_intersect(a, b):
return from_multiset(aa & bb) return from_multiset(aa & bb)
def HMCombine(data, fname=None, type=1, **kwargs): def hm_combine(data, fname=None, type=1, **kwargs):
""" """
Plot heatmap for calibration report Plot heatmap for calibration report
:param data: data to plot :param data: 2D numpy.array of data points to plot
:param fname: file name of output file :param fname: file name of output file
:param type: type of figure :param type: type of figure
:param kwargs: other keywords supported by pyDetLib heatmapPlot :param kwargs: other keywords supported by pyDetLib heatmapPlot
...@@ -334,7 +335,7 @@ def HMCombine(data, fname=None, type=1, **kwargs): ...@@ -334,7 +335,7 @@ def HMCombine(data, fname=None, type=1, **kwargs):
vmin = kwargs.get('vmin', None) vmin = kwargs.get('vmin', None)
vmax = kwargs.get('vmax', None) vmax = kwargs.get('vmax', None)
h_frame = (1 - pad_b - pad_t) h_frame = 1 - pad_b - pad_t
w_frame = 1 - pad_l - pad_r w_frame = 1 - pad_l - pad_r
if type == 2: if type == 2:
...@@ -347,7 +348,8 @@ def HMCombine(data, fname=None, type=1, **kwargs): ...@@ -347,7 +348,8 @@ def HMCombine(data, fname=None, type=1, **kwargs):
mean = np.nanmean(data[y]) mean = np.nanmean(data[y])
if vmin and vmax: if vmin and vmax:
mean = np.nanmean( mean = np.nanmean(
data[y, np.where((data[y] > vmin) & (data[y] < vmax))]) data[y, np.where((data[y] > vmin) & (data[y] < vmax))]
)
ax = plt.axes( ax = plt.axes(
[pad_l, pad_b + h_frame * y + mar, w_frame, h_frame - 2 * mar], [pad_l, pad_b + h_frame * y + mar, w_frame, h_frame - 2 * mar],
...@@ -356,17 +358,16 @@ def HMCombine(data, fname=None, type=1, **kwargs): ...@@ -356,17 +358,16 @@ def HMCombine(data, fname=None, type=1, **kwargs):
yticklabels=['{:4.2f}'.format(mean)], yticklabels=['{:4.2f}'.format(mean)],
xticks=[], xticks=[],
xticklabels=[] xticklabels=[]
) )
ax.tick_params(axis='y', which='major', pad=25) ax.tick_params(axis='y', which='major', pad=25)
if vmin and vmax: if vmin and vmax:
ax.set_ylim([vmin, vmax]) ax.set_ylim(vmin, vmax)
d = [{'x': np.arange(data.shape[1]), d = [{'x': np.arange(data.shape[1]),
'y': data[y], 'y': data[y],
'drawstyle': 'steps-mid', 'drawstyle': 'steps-mid',
'linestyle': 'dotted', 'linestyle': 'dotted',
# 'marker': '^',
'linewidth': 2, 'linewidth': 2,
'color': 'white', 'color': 'white',
...@@ -381,15 +382,16 @@ def HMCombine(data, fname=None, type=1, **kwargs): ...@@ -381,15 +382,16 @@ def HMCombine(data, fname=None, type=1, **kwargs):
] ]
y_lab = '' y_lab = ''
if y == data.shape[0] / 2: # Locate label in a middle of the axis
if y == int(data.shape[0] / 2):
y_lab = 'Average over time' y_lab = 'Average over time'
_ = xana.simplePlot(d, # aspect=1.6, _ = xana.simplePlot(d,
x_label="", x_label="",
y_label=y_lab, y_label=y_lab,
use_axis=ax, use_axis=ax,
y_log=False) y_log=False)
elif type == 1: if type == 1:
ax2 = plt.axes([pad_l, pad_b, w_frame, h_frame / 16.], frame_on=False, ax2 = plt.axes([pad_l, pad_b, w_frame, h_frame / 16.], frame_on=False,
yticks=range(3), yticklabels=[32, 16, 0], yticks=range(3), yticklabels=[32, 16, 0],
...@@ -399,8 +401,6 @@ def HMCombine(data, fname=None, type=1, **kwargs): ...@@ -399,8 +401,6 @@ def HMCombine(data, fname=None, type=1, **kwargs):
ax2.set_ylabel('Memory cell ID', color='r') ax2.set_ylabel('Memory cell ID', color='r')
ax2.tick_params('y', colors='r', right=True, left=False, ax2.tick_params('y', colors='r', right=True, left=False,
labelright=True, labelleft=False, ) labelright=True, labelleft=False, )
else:
pass
if fname: if fname:
fig.savefig(fname, bbox_inches='tight') fig.savefig(fname, bbox_inches='tight')
notebooks/AGIPD/PlotFromCalDB_AGIPD_NBC.ipynb
+ 59
− 81
View file @ bff91e27
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Statistical analysis of calibration factors# # Statistical analysis of calibration factors#
Author: Mikhail Karnevskiy, Steffen Hauf, Version 0.2 Author: Mikhail Karnevskiy, Steffen Hauf, Version 0.2
Calibration constants for AGIPD1M1 detector from the data base with injection time between start_date and end_date are considered. Calibration constants for AGIPD1M1 detector from the data base with injection time between start_date and end_date are considered.
To be visualized, calibration constants are averaged per ASICs. Plots shows calibration constant over time for each constant and for each module. Summary plots overall modules are created. To be visualized, calibration constants are averaged per ASICs. Plots shows calibration constant over time for each constant and for each module. Summary plots overall modules are created.
In additional gain-slopes flat-field and pulse-capacitor are combined to relative-gain constant and presented as well. Noise in electron units is derived using gain factors and presented. In additional gain-slopes flat-field and pulse-capacitor are combined to relative-gain constant and presented as well. Noise in electron units is derived using gain factors and presented.
Values shown in plots are saved in h5 files. Values shown in plots are saved in h5 files.
All presented values corresponds to high and medium gain stages. All presented values corresponds to high and medium gain stages.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
cluster_profile = "noDB" # The ipcluster profile to use cluster_profile = "noDB" # The ipcluster profile to use
start_date = "2018-01-30" # Date to start investigation interval from start_date = "2018-01-30" # Date to start investigation interval from
end_date = "2018-12-12" # Date to end investigation interval at, can be "now" end_date = "2018-12-12" # Date to end investigation interval at, can be "now"
constants = ["Noise", "SlopesFF", "SlopesPC", "Offset"] # Constants to plot constants = ["Noise", "SlopesFF", "SlopesPC", "Offset"] # Constants to plot
modules = [0] # Modules, range allowed modules = [0] # Modules, range allowed
bias_voltages = [300, 500] # Bias voltage bias_voltages = [300, 500] # Bias voltage
mem_cells = [128, 176] # Number of used memory cells. Typically: 4,32,64,128,176. mem_cells = [128, 176] # Number of used memory cells. Typically: 4,32,64,128,176.
photon_energy = 9.2 # Photon energy of the beam photon_energy = 9.2 # Photon energy of the beam
out_folder = "/gpfs/exfel/data/scratch/karnem/testAGIPD16_14/" # Output folder, required out_folder = "/gpfs/exfel/data/scratch/karnem/testAGIPD16_10/" # Output folder, required
use_existing = "/gpfs/exfel/data/scratch/karnem/testAGIPD16_9/" # Input folder use_existing = "" # Input folder
cal_db_timeout = 1200000 # timeout on caldb requests", cal_db_timeout = 1200000 # timeout on caldb requests",
detectorDB = "AGIPD1M1" # detector entry in the DB to investigate detectorDB = "AGIPD1M1" # detector entry in the DB to investigate
dclass = "AGIPD" # Detector class dclass = "AGIPD" # Detector class
cal_db_interface = "tcp://max-exfl016:8015#8025" # the database interface to use cal_db_interface = "tcp://max-exfl016:8015#8025" # the database interface to use
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import os
import numpy as np
import matplotlib
##matplotlib.use("agg")
import matplotlib.pyplot as plt
% matplotlib inline
import sys
import datetime import datetime
import dateutil.parser import dateutil.parser
import numpy as np
import os
import sys
import warnings import warnings
warnings.filterwarnings('ignore') warnings.filterwarnings('ignore')
import matplotlib
% matplotlib inline
from iCalibrationDB import Constants, Conditions, Detectors from iCalibrationDB import Constants, Conditions, Detectors
from cal_tools.tools import get_from_db from cal_tools.tools import get_from_db
from cal_tools.ana_tools import * from cal_tools.ana_tools import *
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Prepare variables # Prepare variables
interval = 1 # interval for evaluation in days interval = 1 # interval for evaluation in days
if modules[0] == -1: nMem = 176 # Number of mem Cells to store
modules = list(range(16)) spShape = (64,64) # Shape of superpixel
adu_to_photon = 8000 / 3.6 / 67.# ADU to photon conversion factor
in_mod_names = [] in_mod_names = []
for mod in modules: for mod in modules:
in_mod_names.append("Q{}M{}".format(mod // 4 + 1, mod % 4 + 1)) in_mod_names.append("Q{}M{}".format(mod // 4 + 1, mod % 4 + 1))
constantsDark = {"SlopesFF": 'BadPixelsFF', constantsDark = {"SlopesFF": 'BadPixelsFF',
'SlopesPC': 'BadPixelsPC', 'SlopesPC': 'BadPixelsPC',
'Noise': 'BadPixelsDark', 'Noise': 'BadPixelsDark',
'Offset': 'BadPixelsDark'} 'Offset': 'BadPixelsDark'}
print('Bad pixels data: ', constantsDark) print('Bad pixels data: ', constantsDark)
# Define parameters in order to perform loop over time stamps # Define parameters in order to perform loop over time stamps
start = datetime.datetime.now() if start_date.upper() == "NOW" else dateutil.parser.parse( start = datetime.datetime.now() if start_date.upper() == "NOW" else dateutil.parser.parse(
start_date) start_date)
end = datetime.datetime.now() if end_date.upper() == "NOW" else dateutil.parser.parse( end = datetime.datetime.now() if end_date.upper() == "NOW" else dateutil.parser.parse(
end_date) end_date)
step = datetime.timedelta(hours=interval) step = datetime.timedelta(hours=interval)
dt = end dt = end
# Create output folder # Create output folder
os.makedirs(out_folder, exist_ok=True) os.makedirs(out_folder, exist_ok=True)
# Get getector conditions # Get getector conditions
det = getattr(Detectors, detectorDB) det = getattr(Detectors, detectorDB)
dconstants = getattr(Constants, dclass) dconstants = getattr(Constants, dclass)
if "#" in cal_db_interface: if "#" in cal_db_interface:
prot, serv, ran = cal_db_interface.split(":") prot, serv, ran = cal_db_interface.split(":")
r1, r2 = ran.split("#") r1, r2 = ran.split("#")
cal_db_interface = ":".join( cal_db_interface = ":".join(
[prot, serv, str(np.random.randint(int(r1), int(r2)))]) [prot, serv, str(np.random.randint(int(r1), int(r2)))])
print('CalDB Interface {}'.format(cal_db_interface)) print('CalDB Interface {}'.format(cal_db_interface))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
```
%% Cell type:code id: tags:
``` python
import copy import copy
def prepare_to_store(a, nMem): def prepare_to_store(a, nMem):
shape = list(a.shape[:2])+[nMem,2] shape = list(a.shape[:2])+[nMem, 2]
b = np.full(shape, np.nan) b = np.full(shape, np.nan)
b[:, :, :a.shape[2]] = a[:, :, :, :2] b[:, :, :a.shape[2]] = a[:, :, :, :2]
return b return b
def get_rebined(a, rebin): def get_rebined(a, rebin):
return a.reshape( return a.reshape(
int(a.shape[0] / rebin[0]), int(a.shape[0] / rebin[0]),
rebin[0], rebin[0],
int(a.shape[1] / rebin[1]), int(a.shape[1] / rebin[1]),
rebin[1], rebin[1],
a.shape[2], a.shape[2],
a.shape[3]) a.shape[3])
def modify_const(const, data, isBP = False):
if const in ['SlopesFF']:
if (len(data.shape) == 4):
data = data[:, :, :, 0][..., None]
else:
data = data[..., None]
nMem = 176 # Number of mem Cells to store if not isBP:
spShape = (64,64) # Shape of superpixel if data.shape[0] != 128:
data = data.swapaxes(0, 2).swapaxes(1, 3).swapaxes(2, 3)
# Copy slope medium to be saved later
if const in ['SlopesPC']:
data[:, :, :, 1] = data[:, :, :, 3]
else:
if const in ['SlopesPC']:
if len(data.shape) == 3:
data = data[:, :, :, None].repeat(10, axis=3)
if data.shape[0] != 128:
data = data.swapaxes(0, 1).swapaxes(1, 2)
if len(data.shape) < 4:
print(data.shape, "Unexpected shape !!!!!!!!!")
return data
ret_constants = {} ret_constants = {}
if use_existing != '': if use_existing != '':
mem_cells = [] mem_cells = []
# Loop over max_memory cells # Loop over max_memory cells
for the_mem_cells in mem_cells: for the_mem_cells in mem_cells:
# Loop over bias voltages # Loop over bias voltages
for bias_voltage in bias_voltages: for bias_voltage in bias_voltages:
# Loop over constants # Loop over constants
for c, const in enumerate(constants): for c, const in enumerate(constants):
if not const in ret_constants: if not const in ret_constants:
ret_constants[const] = {} ret_constants[const] = {}
# Loop over modules # Loop over modules
for module in modules: for module in modules:
dt = end dt = end
qm = "Q{}M{}".format(module // 4 + 1, module % 4 + 1) qm = "Q{}M{}".format(module // 4 + 1, module % 4 + 1)
if not qm in ret_constants[const]: if not qm in ret_constants[const]:
ret_constants[const][qm] = [] ret_constants[const][qm] = []
# Loop over time stamps # Loop over time stamps
while dt > start: while dt > start:
creation_time = dt creation_time = dt
if (const in ["Offset", "Noise", if (const in ["Offset", "Noise",
"SlopesPC"] or "DARK" in const.upper()): "SlopesPC"] or "DARK" in const.upper()):
dcond = Conditions.Dark dcond = Conditions.Dark
mcond = getattr(dcond, dclass)( mcond = getattr(dcond, dclass)(
memory_cells=the_mem_cells, memory_cells=the_mem_cells,
bias_voltage=bias_voltage) bias_voltage=bias_voltage)
else: else:
dcond = Conditions.Illuminated dcond = Conditions.Illuminated
mcond = getattr(dcond, dclass)( mcond = getattr(dcond, dclass)(
memory_cells=the_mem_cells, memory_cells=the_mem_cells,
bias_voltage=bias_voltage, bias_voltage=bias_voltage,
photon_energy=photon_energy) photon_energy=photon_energy)
print('Request: ', const, qm, the_mem_cells, bias_voltage, print('Request: ', const, qm, the_mem_cells, bias_voltage,
creation_time) creation_time)
cdata, ctime = get_from_db(getattr(det, qm), cdata, ctime = get_from_db(getattr(det, qm),
getattr(dconstants, const)(), getattr(dconstants, const)(),
copy.deepcopy(mcond), copy.deepcopy(mcond),
None, None,
cal_db_interface, cal_db_interface,
creation_time=creation_time, creation_time=creation_time,
verbosity=0, verbosity=0,
timeout=cal_db_timeout, timeout=cal_db_timeout,
meta_only=True ) meta_only=True )
if ctime is None or cdata is None: if ctime is None or cdata is None:
print('Time or Data is None') print('Time or Data is None')
break break
ctime = ctime.calibration_constant_version.begin_at ctime = ctime.calibration_constant_version.begin_at
print("Found constant {}: begin at {}".format(const, print("Found constant {}: begin at {}".format(const,
cdata is not None), cdata is not None),
ctime) ctime)
print(cdata.shape) print(cdata.shape)
cdata = modify_const(const, cdata)
if const in ['SlopesFF']:
if (len(cdata.shape) == 4):
cdata = cdata[:, :, :, 0][..., None]
else:
cdata = cdata[..., None]
if len(cdata.shape) < 4:
print(cdata.shape, "Unexpected shape !!!!!!!!!")
if cdata.shape[0] != 128:
cdata = cdata.swapaxes(0, 2).swapaxes(1, 3).swapaxes(2, 3)
# Copy slope medium to be saved later
if const in ['SlopesPC']:
cdata[:, :, :, 1] = cdata[:, :, :, 3]
print(cdata.shape) print(cdata.shape)
# Request bad pixel mask # Request bad pixel mask
if const in constantsDark: if const in constantsDark:
#for par in mcond.parameters:
# print (par.name, par.value)
print('RequestBP ', print('RequestBP ',
(creation_time + step + step + step + step)) (creation_time + step + step + step + step))
cdataBP, ctimeBP = get_from_db(getattr(det, qm), cdataBP, ctimeBP = get_from_db(getattr(det, qm),
getattr(dconstants, getattr(dconstants,
constantsDark[const])(), constantsDark[const])(),
copy.deepcopy(mcond), copy.deepcopy(mcond),
None, None,
cal_db_interface, cal_db_interface,
creation_time=( creation_time=(
creation_time + step + step + step + step), creation_time + step + step + step + step),
verbosity=0, verbosity=0,
timeout=cal_db_timeout, timeout=cal_db_timeout,
meta_only=True) meta_only=True)
ctimeBP = ctimeBP.calibration_constant_version.begin_at ctimeBP = ctimeBP.calibration_constant_version.begin_at
if cdataBP is None: if cdataBP is None:
print("Bad pixels are not found !!!!!!!!!") print("Bad pixels are not found !!!!!!!!!")
ctime = ctime.replace(tzinfo=None) ctime = ctime.replace(tzinfo=None)
dt = ctime - step dt = ctime - step
continue continue
print("Found BP {}".format(ctimeBP)) print("Found BP {}".format(ctimeBP))
print(cdataBP.shape) print(cdataBP.shape)
cdataBP = modify_const(const, cdataBP, True)
if const in ['SlopesFF']: print(cdataBP.shape)
if len(cdataBP.shape) == 4:
cdataBP = cdataBP[:, :, :, 0][..., None]
else:
cdataBP = cdataBP[..., None]
if const in ['SlopesPC']:
if len(cdataBP.shape) == 3:
cdataBP = cdataBP[:, :, :, None].repeat(10,
axis=3)
if cdataBP.shape[0] != 128:
cdataBP = cdataBP.swapaxes(0, 1).swapaxes(1, 2)
if (len(cdataBP.shape) < 4):
print(cdataBP.shape, "Unexpected shape !!!!!!!!!")
if cdataBP.shape != cdata.shape: if cdataBP.shape != cdata.shape:
print(cdataBP.shape, 'Wrong BP shape!!!') print('Wrong BP shape!!!')
ctime = ctime.replace(tzinfo=None) ctime = ctime.replace(tzinfo=None)
dt = ctime - step dt = ctime - step
continue continue
# Apply bad pixel mask # Apply bad pixel mask
cdataABP = np.copy(cdata) cdataABP = np.copy(cdata)
cdataABP[cdataBP > 0] = np.nan cdataABP[cdataBP > 0] = np.nan
# Create superpixels for constants with BP applied # Create superpixels for constants with BP applied
cdataABP = get_rebined(cdataABP, spShape) cdataABP = get_rebined(cdataABP, spShape)
toStoreBP = prepare_to_store(np.nanmean(cdataABP, axis=(1, 3)), nMem) toStoreBP = prepare_to_store(np.nanmean(cdataABP, axis=(1, 3)), nMem)
toStoreBPStd = prepare_to_store(np.nanstd(cdataABP, axis=(1, 3)), nMem) toStoreBPStd = prepare_to_store(np.nanstd(cdataABP, axis=(1, 3)), nMem)
# Prepare number of bad pixels per superpixels # Prepare number of bad pixels per superpixels
cdataBP = get_rebined(cdataBP, spShape) cdataBP = get_rebined(cdataBP, spShape)
cdataNBP = prepare_to_store(np.nansum(cdataBP > 0, axis=(1, 3)), nMem) cdataNBP = prepare_to_store(np.nansum(cdataBP > 0, axis=(1, 3)), nMem)
else: else:
toStoreExtBP = 0 toStoreExtBP = 0
cdataBPExt = 0 cdataBPExt = 0
# Create superpixels for constants without BP applied # Create superpixels for constants without BP applied
cdata = get_rebined(cdata, spShape) cdata = get_rebined(cdata, spShape)
toStoreStd = prepare_to_store(np.nanstd(cdata, axis=(1, 3)), nMem) toStoreStd = prepare_to_store(np.nanstd(cdata, axis=(1, 3)), nMem)
toStore = prepare_to_store(np.nanmean(cdata, axis=(1, 3)), nMem) toStore = prepare_to_store(np.nanmean(cdata, axis=(1, 3)), nMem)
# Convert parameters to dict # Convert parameters to dict
dpar = {} dpar = {}
for par in mcond.parameters: for par in mcond.parameters:
dpar[par.name] = par.value dpar[par.name] = par.value
print("Store values in dict", const, qm, ctime) print("Store values in dict", const, qm, ctime)
ret_constants[const][qm].append({'ctime': ctime, ret_constants[const][qm].append({'ctime': ctime,
'nBP': cdataNBP, 'nBP': cdataNBP,
'dataBP': toStoreBP, 'dataBP': toStoreBP,
'dataBPStd': toStoreBPStd, 'dataBPStd': toStoreBPStd,
'data': toStore, 'data': toStore,
'dataStd': toStoreStd, 'dataStd': toStoreStd,
'mdata': dpar}) 'mdata': dpar})
ctime = ctime.replace(tzinfo=None) ctime = ctime.replace(tzinfo=None)
dt = ctime - step dt = ctime - step
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if use_existing == "": if use_existing == "":
print('Save data to /CalDBAna_{}_{}.h5'.format(dclass, modules[0])) print('Save data to /CalDBAna_{}_{}.h5'.format(dclass, modules[0]))
save_dict_to_hdf5(ret_constants, save_dict_to_hdf5(ret_constants,
'{}/CalDBAna_{}_{}.h5'.format(out_folder, dclass, modules[0])) '{}/CalDBAna_{}_{}.h5'.format(out_folder, dclass, modules[0]))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if use_existing == "": if use_existing == "":
fpath = '{}/CalDBAna_{}_*.h5'.format(out_folder, dclass) fpath = '{}/CalDBAna_{}_*.h5'.format(out_folder, dclass)
else: else:
fpath = '{}/CalDBAna_{}_*.h5'.format(use_existing, dclass) fpath = '{}/CalDBAna_{}_*.h5'.format(use_existing, dclass)
print('Load data from {}'.format(fpath)) print('Load data from {}'.format(fpath))
ret_constants = load_data_from_hdf5(fpath) ret_constants = load_data_from_hdf5(fpath)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
print ('Combine Gain FF and PC and Noise in [e-]') # Combine FF and PC data to calculate Gain
# Estimate Noise in units of electrons
print ('Calculate Gain and Noise in electron units')
ret_constants["Gain"] = {} ret_constants["Gain"] = {}
ret_constants["Noise-e"] = {} ret_constants["Noise-e"] = {}
for module in list(range(16)): for module in list(range(16)):
if ("SlopesFF" not in ret_constants or if ("SlopesFF" not in ret_constants or
"SlopesPC" not in ret_constants): "SlopesPC" not in ret_constants):
break break
qm = "Q{}M{}".format(module // 4 + 1, module % 4 + 1) qm = "Q{}M{}".format(module // 4 + 1, module % 4 + 1)
print(qm) print(qm)
if (qm not in ret_constants["SlopesFF"] or if (qm not in ret_constants["SlopesFF"] or
qm not in ret_constants["SlopesPC"]): qm not in ret_constants["SlopesPC"]):
continue continue
ret_constants["Gain"][qm] = {} ret_constants["Gain"][qm] = {}
dataFF = ret_constants["SlopesFF"][qm] dataFF = ret_constants["SlopesFF"][qm]
dataPC = ret_constants["SlopesPC"][qm] dataPC = ret_constants["SlopesPC"][qm]
if (len(dataFF) == 0 or len(dataPC) == 0): if (len(dataFF) == 0 or len(dataPC) == 0):
continue continue
ctimesFF = np.array(dataFF["ctime"]) ctimesFF = np.array(dataFF["ctime"])
ctimesPC = np.array(dataPC["ctime"]) ctimesPC = np.array(dataPC["ctime"])
ctime, icomb = combine_constants(ctimesFF, ctimesPC) ctime, icomb = combine_constants(ctimesFF, ctimesPC)
cdataPC_vs_time = np.array(dataPC["data"])[..., 0] cdataPC_vs_time = np.array(dataPC["data"])[..., 0]
cdataFF_vs_time = np.array(dataFF["data"])[..., 0] cdataFF_vs_time = np.array(dataFF["data"])[..., 0]
cdataFF_vs_time = np.nanmedian(cdataFF_vs_time, axis=3)[..., None] cdataFF_vs_time = np.nanmedian(cdataFF_vs_time, axis=3)[..., None]
cdataFF_vs_time /= np.nanmedian(cdataFF_vs_time, axis=(1, 2, 3))[:, None, cdataFF_vs_time /= np.nanmedian(cdataFF_vs_time, axis=(1, 2, 3))[:, None,
None, None] None, None]
cdataPC_vs_time /= np.nanmedian(cdataPC_vs_time, axis=(1, 2, 3))[:, None, cdataPC_vs_time /= np.nanmedian(cdataPC_vs_time, axis=(1, 2, 3))[:, None,
None, None] None, None]
gain_vs_time = [] gain_vs_time = []
for iFF, iPC in icomb: for iFF, iPC in icomb:
gain_vs_time.append(cdataFF_vs_time[iFF] * cdataPC_vs_time[iPC]) gain_vs_time.append(cdataFF_vs_time[iFF] * cdataPC_vs_time[iPC])
print(np.array(gain_vs_time).shape) print(np.array(gain_vs_time).shape)
ctime_ts = [t.timestamp() for t in ctime] ctime_ts = [t.timestamp() for t in ctime]
ret_constants["Gain"][qm]["ctime"] = ctime ret_constants["Gain"][qm]["ctime"] = ctime
ret_constants["Gain"][qm]["data"] = np.array(gain_vs_time) ret_constants["Gain"][qm]["data"] = np.array(gain_vs_time)
# Fill missing data for compatibility with plotting code
ret_constants["Gain"][qm]["dataBP"] = np.array(gain_vs_time) ret_constants["Gain"][qm]["dataBP"] = np.array(gain_vs_time)
ret_constants["Gain"][qm]["nBP"] = np.array(gain_vs_time) ret_constants["Gain"][qm]["nBP"] = np.array(gain_vs_time)
if "Noise" not in ret_constants: if "Noise" not in ret_constants:
continue continue
if qm not in ret_constants["Noise"]: if qm not in ret_constants["Noise"]:
continue continue
dataN = ret_constants["Noise"][qm] dataN = ret_constants["Noise"][qm]
if len(dataN) == 0: if len(dataN) == 0:
continue continue
ret_constants["Noise-e"][qm] = {} ret_constants["Noise-e"][qm] = {}
ctimesG = np.array(ctime) ctimesG = np.array(ctime)
ctimesN = np.array(dataN["ctime"]) ctimesN = np.array(dataN["ctime"])
ctime, icomb = combine_constants(ctimesG, ctimesN) ctime, icomb = combine_constants(ctimesG, ctimesN)
cdataG_vs_time = np.array(gain_vs_time) cdataG_vs_time = np.array(gain_vs_time)
cdataN_vs_time = np.array(dataN["data"])[..., 0] cdataN_vs_time = np.array(dataN["data"])[..., 0]
data_vs_time = [] data_vs_time = []
for iG, iN in icomb: for iG, iN in icomb:
data_vs_time.append( data_vs_time.append(
cdataN_vs_time[iN] * (8000 / 3.6 / 67.) / cdataG_vs_time[iG]) cdataN_vs_time[iN] * adu_to_photon / cdataG_vs_time[iG])
print(np.array(gain_vs_time).shape) print(np.array(gain_vs_time).shape)
ctime_ts = [t.timestamp() for t in ctime] ctime_ts = [t.timestamp() for t in ctime]
ret_constants["Noise-e"][qm]["ctime"] = ctime ret_constants["Noise-e"][qm]["ctime"] = ctime
ret_constants["Noise-e"][qm]["data"] = np.array(data_vs_time) ret_constants["Noise-e"][qm]["data"] = np.array(data_vs_time)
# Fill missing data for compatibility with plotting code
ret_constants["Noise-e"][qm]["dataBP"] = np.array(data_vs_time) ret_constants["Noise-e"][qm]["dataBP"] = np.array(data_vs_time)
ret_constants["Noise-e"][qm]["nBP"] = np.array(data_vs_time) ret_constants["Noise-e"][qm]["nBP"] = np.array(data_vs_time)
save_dict_to_hdf5({k:v for k,v in ret_constants.items() if k in ['Gain', 'Noise-e']}, save_dict_to_hdf5({k:v for k,v in ret_constants.items() if k in ['Gain', 'Noise-e']},
'{}/CalDBAna_{}_Gain.h5'.format(out_folder, dclass)) '{}/CalDBAna_{}_Gain.h5'.format(out_folder, dclass))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Parameters for plotting # Parameters for plotting
# Define range for plotting # Define range for plotting
rangevals = { rangevals = {
"Offset": [[800., 1500], [600, 900]], "Offset": [[4000., 5500], [6500, 8500]],
"Noise": [[2.0, 16], [1.0, 7.0]], "Noise": [[2.5, 15], [7.5, 17.0]],
"Gain": [[20, 30], [20, 30]], "Gain": [[0.8, 1.2], [0.8, 1.2]],
"Noise-e": [[100., 600.], [100., 600.]], "Noise-e": [[85., 500.], [85., 500.]],
"SlopesCI": [[0.95, 1.05], [0.0, 0.5]], "SlopesPC": [[22.0, 27.0], [-0.5, 1.5]],
"SlopesFF": [[0.8, 1.2], [0.8, 1.2]] "SlopesFF": [[0.8, 1.2], [0.6, 1.2]]
} }
nMemToShow = 32 nMemToShow = 32
keys = { keys = {
'Mean': ['data', '', 'Mean over pixels'], 'Mean': ['data', '', 'Mean over pixels'],
'std': ['dataStd', '', '$\sigma$ over pixels'], 'std': ['dataStd', '', '$\sigma$ over pixels'],
'MeanBP': ['dataBP', 'Good pixels only', 'Mean over pixels'], 'MeanBP': ['dataBP', 'Good pixels only', 'Mean over pixels'],
'NBP': ['nBP', 'Fraction of BP', 'Fraction of BP'], 'NBP': ['nBP', 'Fraction of BP', 'Fraction of BP'],
'stdBP': ['dataBPStd', 'Good pixels only', '$\sigma$ over pixels'], 'stdBP': ['dataBPStd', 'Good pixels only', '$\sigma$ over pixels'],
'stdASIC': ['', '', '$\sigma$ over ASICs'], 'stdASIC': ['', '', '$\sigma$ over ASICs'],
'stdCell': ['', '', '$\sigma$ over Cells'], 'stdCell': ['', '', '$\sigma$ over Cells'],
} }
gain_name = ['High', 'Medium', 'Low'] gain_name = ['High', 'Medium', 'Low']
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
print('Plot calibration constants') print('Plot calibration constants')
# loop over constat type # loop over constat type
for const, modules in ret_constants.items(): for const, modules in ret_constants.items():
# Loop over gain # Loop over gain
for gain in range(2): for gain in range(2):
print('Const: {}, gain {}'.format(const, gain)) print('Const: {}, gain {}'.format(const, gain))
if const in ["Gain", "Noise-e"] and gain == 1: if const in ["Gain", "Noise-e"] and gain == 1:
continue continue
else: else:
pass pass
# loop over modules # loop over modules
mod_data = {} mod_data = {}
mod_data['stdASIC'] = [] mod_data['stdASIC'] = []
mod_data['stdCell'] = [] mod_data['stdCell'] = []
mod_names = [] mod_names = []
mod_times = [] mod_times = []
# Loop over modules # Loop over modules
for mod, data in modules.items(): for mod, data in modules.items():
if mod not in in_mod_names: if mod not in in_mod_names:
continue continue
print(mod) print(mod)
ctimes = np.array(data["ctime"]) ctimes = np.array(data["ctime"])
ctimes_ticks = [x.strftime('%m-%d') for x in ctimes] ctimes_ticks = [x.strftime('%m-%d') for x in ctimes]
if ("mdata" in data): if ("mdata" in data):
cmdata = np.array(data["mdata"]) cmdata = np.array(data["mdata"])
for i, tick in enumerate(ctimes_ticks): for i, tick in enumerate(ctimes_ticks):
ctimes_ticks[i] = ctimes_ticks[i] + \ ctimes_ticks[i] = ctimes_ticks[i] + \
', V={:1.0f}'.format(cmdata[i]['Sensor Bias Voltage']) + \ ', V={:1.0f}'.format(cmdata[i]['Sensor Bias Voltage']) + \
', M={:1.0f}'.format( ', M={:1.0f}'.format(
cmdata[i]['Memory cells']) cmdata[i]['Memory cells'])
sort_ind = np.argsort(ctimes_ticks) sort_ind = np.argsort(ctimes_ticks)
ctimes_ticks = list(np.array(ctimes_ticks)[sort_ind]) ctimes_ticks = list(np.array(ctimes_ticks)[sort_ind])
# Create sorted by data dataset # Create sorted by data dataset
rdata = {} rdata = {}
for key, item in keys.items(): for key, item in keys.items():
if item[0] in data: if item[0] in data:
rdata[key] = np.array(data[item[0]])[sort_ind] rdata[key] = np.array(data[item[0]])[sort_ind]
# print(rdata['Mean'].shape)
nTimes = rdata['Mean'].shape[0] nTimes = rdata['Mean'].shape[0]
nPixels = rdata['Mean'].shape[1] * rdata['Mean'].shape[2] nPixels = rdata['Mean'].shape[1] * rdata['Mean'].shape[2]
nBins = nMemToShow * nPixels nBins = nMemToShow * nPixels
# Select gain # Select gain
if const not in ["Gain", "Noise-e"]: if const not in ["Gain", "Noise-e"]:
for key in rdata: for key in rdata:
rdata[key] = rdata[key][..., gain] rdata[key] = rdata[key][..., gain]
# Avoid to low values # Avoid to low values
if const in ["Noise", "Offset", "Noise-e"]: if const in ["Noise", "Offset", "Noise-e"]:
rdata['Mean'][rdata['Mean'] < 0.1] = np.nan rdata['Mean'][rdata['Mean'] < 0.1] = np.nan
if 'MeanBP' in rdata: if 'MeanBP' in rdata:
rdata['MeanBP'][rdata['MeanBP'] < 0.1] = np.nan rdata['MeanBP'][rdata['MeanBP'] < 0.1] = np.nan
if 'NBP' in rdata: if 'NBP' in rdata:
rdata["NBP"][rdata["NBP"] == 4096] = np.nan rdata["NBP"][rdata["NBP"] == 4096] = np.nan
rdata["NBP"] = rdata["NBP"] / (64 * 64) * 100 rdata["NBP"] = rdata["NBP"] / (64 * 64) * 100
# Reshape: ASICs over cells for plotting # Reshape: ASICs over cells for plotting
pdata = {} pdata = {}
for key in rdata: for key in rdata:
pdata[key] = rdata[key][:, :, :, :nMemToShow].reshape( pdata[key] = rdata[key][:, :, :, :nMemToShow].reshape(
nTimes, nBins).swapaxes(0, 1) nTimes, nBins).swapaxes(0, 1)
# Summary over ASICs # Summary over ASICs
adata = {} adata = {}
for key in rdata: for key in rdata:
adata[key] = np.nanmean(rdata[key], axis=(1, 2)).swapaxes(0, 1) adata[key] = np.nanmean(rdata[key], axis=(1, 2)).swapaxes(0, 1)
# Summary information over modules # Summary information over modules
for key in pdata: for key in pdata:
if key not in mod_data: if key not in mod_data:
mod_data[key] = [] mod_data[key] = []
mod_data[key].append(np.nanmean(pdata[key], axis=0)) mod_data[key].append(np.nanmean(pdata[key], axis=0))
mod_data['stdASIC'].append(np.nanstd( mod_data['stdASIC'].append(np.nanstd(
np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=(1, 2)), axis=1)) np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=(1, 2)), axis=1))
mod_data['stdCell'].append(np.nanstd( mod_data['stdCell'].append(np.nanstd(
np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=3), axis=(1, 2))) np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=3), axis=(1, 2)))
mod_names.append(mod) mod_names.append(mod)
mod_times.append(ctimes_ticks) mod_times.append(ctimes_ticks)
# Plotting # Plotting
for key in pdata: for key in pdata:
vmin = None vmin = None
vmax = None vmax = None
if const in rangevals and key in ['Mean', 'MeanBP']: if const in rangevals and key in ['Mean', 'MeanBP']:
vmin = rangevals[const][gain][0] vmin = rangevals[const][gain][0]
vmax = rangevals[const][gain][1] vmax = rangevals[const][gain][1]
if key == 'NBP': if key == 'NBP':
unit = '[%]' unit = '[%]'
else: else:
unit = '[ADU]' unit = '[ADU]'
if const == 'Noise-e': if const == 'Noise-e':
unit = '[$e^-$]' unit = '[$e^-$]'
title = '{}, module {}, {} gain, {}'.format( title = '{}, module {}, {} gain, {}'.format(
const, mod, gain_name[gain], keys[key][1]) const, mod, gain_name[gain], keys[key][1])
cb_label = '{}, {} {}'.format(const, keys[key][2], unit) cb_label = '{}, {} {}'.format(const, keys[key][2], unit)
HMCombine(pdata[key][::-1], HMCombine(pdata[key][::-1],
x_label='Creation Time', y_label='ASIC ID', x_label='Creation Time', y_label='ASIC ID',
x_ticklabels=ctimes_ticks, x_ticklabels=ctimes_ticks,
x_ticks=np.arange(len(ctimes_ticks))+0.3, x_ticks=np.arange(len(ctimes_ticks))+0.3,
title=title, cb_label=cb_label, title=title, cb_label=cb_label,
vmin=vmin, vmax=vmax, vmin=vmin, vmax=vmax,
fname='{}/{}_{}_g{}_ASIC_{}.png'.format( fname='{}/{}_{}_g{}_ASIC_{}.png'.format(
out_folder, const, mod, gain, key), out_folder, const, mod, gain, key),
y_ticks=np.arange(nBins, step=nMemToShow)+16, y_ticks=np.arange(nBins, step=nMemToShow)+16,
y_ticklabels=np.arange(nPixels)[::-1]+1, y_ticklabels=np.arange(nPixels)[::-1]+1,
pad=[0.125, 0.125, 0.12, 0.185]) pad=[0.125, 0.125, 0.12, 0.185])
HMCombine(adata[key], type=0, HMCombine(adata[key], type=0,
x_label='Creation Time', y_label='Memory cell ID', x_label='Creation Time', y_label='Memory cell ID',
x_ticklabels=ctimes_ticks, x_ticklabels=ctimes_ticks,
x_ticks=np.arange(len(ctimes_ticks))+0.3, x_ticks=np.arange(len(ctimes_ticks))+0.3,
title=title, cb_label=cb_label, title=title, cb_label=cb_label,
fname='{}/{}_{}_g{}_MEM_{}.png'.format( fname='{}/{}_{}_g{}_MEM_{}.png'.format(
out_folder, const, mod, gain, key), out_folder, const, mod, gain, key),
vmin=vmin, vmax=vmax) vmin=vmin, vmax=vmax)
# break
#break
#break
``` ```
......
notebooks/AGIPD/PlotFromCalDB_Summary_AGIPD_NBC.ipynb
+ 6
− 12
View file @ bff91e27
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Statistical analysis of calibration factors# # Statistical analysis of calibration factors#
Author: Mikhail Karnevskiy, Steffen Hauf, Version 0.2 Author: Mikhail Karnevskiy, Steffen Hauf, Version 0.2
Plot calibration constants for AGIPD1M1 detector aggregated within detector modules. Input information is taken from folder `use_existing`. Corresponding files are prepared by `PlotFromCalDB` notebook. Plot calibration constants for AGIPD1M1 detector aggregated within detector modules. Input information is taken from folder `use_existing`. Corresponding files are prepared by `PlotFromCalDB` notebook.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
cluster_profile = "noDB" # The ipcluster profile to use cluster_profile = "noDB" # The ipcluster profile to use
out_folder = "/gpfs/exfel/data/scratch/karnem/testAGIPD16_14/" # Output folder, required out_folder = "/gpfs/exfel/data/scratch/karnem/testAGIPD16_14/" # Output folder, required
use_existing = "/gpfs/exfel/data/scratch/karnem/testAGIPD16_21/" # Input folder use_existing = "/gpfs/exfel/data/scratch/karnem/testAGIPD16_21/" # Input folder
dclass = "AGIPD" # Detector class dclass = "AGIPD" # Detector class
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from cal_tools.ana_tools import * import warnings
warnings.filterwarnings('ignore')
import numpy as np import numpy as np
import matplotlib import matplotlib
matplotlib.use("agg")
% matplotlib inline % matplotlib inline
import matplotlib.pyplot as plt
import warnings from cal_tools.ana_tools import *
warnings.filterwarnings('ignore')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if use_existing == '': if use_existing == '':
use_existing = out_folder use_existing = out_folder
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
print('Load data from {}/CalDBAna_{}_*.h5'.format(use_existing, dclass)) print('Load data from {}/CalDBAna_{}_*.h5'.format(use_existing, dclass))
ret_constants = load_data_from_hdf5( ret_constants = load_data_from_hdf5(
'{}/CalDBAna_{}_*.h5'.format(use_existing, dclass)) '{}/CalDBAna_{}_*.h5'.format(use_existing, dclass))
print('Evaluated data from {}/CalDBAna_{}_Gain*.h5'.format(out_folder, dclass)) print('Evaluated data from {}/CalDBAna_{}_Gain*.h5'.format(out_folder, dclass))
ret_constants.update(load_data_from_hdf5( ret_constants.update(load_data_from_hdf5(
'{}/CalDBAna_{}_Gain*.h5'.format(out_folder, dclass))) '{}/CalDBAna_{}_Gain*.h5'.format(out_folder, dclass)))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Parameters for plotting # Parameters for plotting
# Define range for plotting # Define range for plotting
rangevals = { rangevals = {
"Offset": [[4000., 5500], [6500, 8500]], "Offset": [[4000., 5500], [6500, 8500]],
"Noise": [[2.5, 15], [7.5, 17.0]], "Noise": [[2.5, 15], [7.5, 17.0]],
"Gain": [[0.8, 1.2], [0.8, 1.2]], "Gain": [[0.8, 1.2], [0.8, 1.2]],
"Noise-e": [[85., 500.], [85., 500.]], "Noise-e": [[85., 500.], [85., 500.]],
"SlopesCI": [[0.95, 1.05], [0.0, 0.5]], "SlopesCI": [[0.95, 1.05], [0.0, 0.5]],
"SlopesPC": [[22.0, 27.0], [-0.5, 1.5]], "SlopesPC": [[22.0, 27.0], [-0.5, 1.5]],
"SlopesFF": [[0.8, 1.2], [0.6, 1.2]] "SlopesFF": [[0.8, 1.2], [0.6, 1.2]]
} }
nMemToShow = 32 nMemToShow = 32
keys = { keys = {
'Mean': ['data', '', 'Mean over pixels'], 'Mean': ['data', '', 'Mean over pixels'],
'std': ['dataStd', '', '$\sigma$ over pixels'], 'std': ['dataStd', '', '$\sigma$ over pixels'],
'MeanBP': ['dataBP', 'Good pixels only', 'Mean over pixels'], 'MeanBP': ['dataBP', 'Good pixels only', 'Mean over pixels'],
'NBP': ['nBP', 'Fraction of BP', 'Fraction of BP'], 'NBP': ['nBP', 'Fraction of BP', 'Fraction of BP'],
'stdBP': ['dataBPStd', 'Good pixels only', '$\sigma$ over pixels'], 'stdBP': ['dataBPStd', 'Good pixels only', '$\sigma$ over pixels'],
'stdASIC': ['', '', '$\sigma$ over ASICs'], 'stdASIC': ['', '', '$\sigma$ over ASICs'],
'stdCell': ['', '', '$\sigma$ over Cells'], 'stdCell': ['', '', '$\sigma$ over Cells'],
} }
gain_name = ['High', 'Medium', 'Low'] gain_name = ['High', 'Medium', 'Low']
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
print('Plot calibration constants') print('Plot calibration constants')
# loop over constat type # loop over constat type
for const, modules in ret_constants.items(): for const, modules in ret_constants.items():
# Loop over gain # Loop over gain
for gain in range(2): for gain in range(2):
print('Const: {}, gain : {}'.format(const, gain)) print('Const: {}, gain : {}'.format(const, gain))
if const in ["Gain", "Noise-e"] and gain == 1: if const in ["Gain", "Noise-e"] and gain == 1:
continue continue
# loop over modules # loop over modules
mod_data = {} mod_data = {}
mod_data['stdASIC'] = [] mod_data['stdASIC'] = []
mod_data['stdCell'] = [] mod_data['stdCell'] = []
mod_names = [] mod_names = []
mod_times = [] mod_times = []
# Loop over modules # Loop over modules
for mod, data in modules.items(): for mod, data in modules.items():
ctimes = np.array(data["ctime"]) ctimes = np.array(data["ctime"])
ctimes_ticks = [x.strftime('%y-%m-%d') for x in ctimes] ctimes_ticks = [x.strftime('%y-%m-%d') for x in ctimes]
if ("mdata" in data): if ("mdata" in data):
cmdata = np.array(data["mdata"]) cmdata = np.array(data["mdata"])
for i, tick in enumerate(ctimes_ticks): for i, tick in enumerate(ctimes_ticks):
ctimes_ticks[i] = ctimes_ticks[i] + \ ctimes_ticks[i] = ctimes_ticks[i] + \
', V={:1.0f}'.format(cmdata[i]['Sensor Bias Voltage']) + \ ', V={:1.0f}'.format(cmdata[i]['Sensor Bias Voltage']) + \
', M={:1.0f}'.format( ', M={:1.0f}'.format(
cmdata[i]['Memory cells']) cmdata[i]['Memory cells'])
sort_ind = np.argsort(ctimes_ticks) sort_ind = np.argsort(ctimes_ticks)
ctimes_ticks = list(np.array(ctimes_ticks)[sort_ind]) ctimes_ticks = list(np.array(ctimes_ticks)[sort_ind])
# Create sorted by data dataset # Create sorted by data dataset
rdata = {} rdata = {}
for key, item in keys.items(): for key, item in keys.items():
if item[0] in data: if item[0] in data:
rdata[key] = np.array(data[item[0]])[sort_ind] rdata[key] = np.array(data[item[0]])[sort_ind]
nTimes = rdata['Mean'].shape[0] nTimes = rdata['Mean'].shape[0]
nPixels = rdata['Mean'].shape[1] * rdata['Mean'].shape[2] nPixels = rdata['Mean'].shape[1] * rdata['Mean'].shape[2]
nBins = nMemToShow * nPixels nBins = nMemToShow * nPixels
# Select gain # Select gain
if const not in ["Gain", "Noise-e"]: if const not in ["Gain", "Noise-e"]:
for key in rdata: for key in rdata:
rdata[key] = rdata[key][..., gain] rdata[key] = rdata[key][..., gain]
# Avoid to low values # Avoid to low values
if const in ["Noise", "Offset", "Noise-e"]: if const in ["Noise", "Offset", "Noise-e"]:
rdata['Mean'][rdata['Mean'] < 0.1] = np.nan rdata['Mean'][rdata['Mean'] < 0.1] = np.nan
if 'MeanBP' in rdata: if 'MeanBP' in rdata:
rdata['MeanBP'][rdata['MeanBP'] < 0.1] = np.nan rdata['MeanBP'][rdata['MeanBP'] < 0.1] = np.nan
if 'NBP' in rdata: if 'NBP' in rdata:
rdata["NBP"][rdata["NBP"] == 4096] = np.nan rdata["NBP"][rdata["NBP"] == 4096] = np.nan
rdata["NBP"] = rdata["NBP"] / (64 * 64) * 100 rdata["NBP"] = rdata["NBP"] / (64 * 64) * 100
# Reshape: ASICs over cells for plotting # Reshape: ASICs over cells for plotting
pdata = {} pdata = {}
for key in rdata: for key in rdata:
pdata[key] = rdata[key][:, :, :, :nMemToShow].reshape( pdata[key] = rdata[key][:, :, :, :nMemToShow].reshape(
nTimes, nBins).swapaxes(0, 1) nTimes, nBins).swapaxes(0, 1)
# Summary information over modules # Summary information over modules
for key in pdata: for key in pdata:
if key not in mod_data: if key not in mod_data:
mod_data[key] = [] mod_data[key] = []
mod_data[key].append(np.nanmean(pdata[key], axis=0)) mod_data[key].append(np.nanmean(pdata[key], axis=0))
mod_data['stdASIC'].append(np.nanstd( mod_data['stdASIC'].append(np.nanstd(
np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=(1, 2)), axis=1)) np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=(1, 2)), axis=1))
mod_data['stdCell'].append(np.nanstd( mod_data['stdCell'].append(np.nanstd(
np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=3), axis=(1, 2))) np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=3), axis=(1, 2)))
mod_names.append(mod) mod_names.append(mod)
mod_times.append(ctimes_ticks) mod_times.append(ctimes_ticks)
# Incert nans to get array-like list of data # Incert nans to get array-like list of data
uTime = mod_times[0] uTime = mod_times[0]
for tlist in mod_times: for tlist in mod_times:
uTime = sorted(multi_union(uTime, tlist)) uTime = sorted(multi_union(uTime, tlist))
for i, tlist in enumerate(mod_times): for i, tlist in enumerate(mod_times):
for t, time in enumerate(uTime): for t, time in enumerate(uTime):
if t == len(tlist) or time != tlist[t]: if t == len(tlist) or time != tlist[t]:
tlist.insert(t, time) tlist.insert(t, time)
for key in mod_data: for key in mod_data:
mod_data[key][i] = np.insert( mod_data[key][i] = np.insert(
mod_data[key][i], t, np.nan) mod_data[key][i], t, np.nan)
# Plotting # Plotting
nModules = len(mod_names) nModules = len(mod_names)
mod_idx = np.argsort(mod_names) mod_idx = np.argsort(mod_names)
for key in mod_data: for key in mod_data:
vmin = None vmin = None
vmax = None vmax = None
if const in rangevals and key in ['Mean', 'MeanBP']: if const in rangevals and key in ['Mean', 'MeanBP']:
vmin = rangevals[const][gain][0] vmin = rangevals[const][gain][0]
vmax = rangevals[const][gain][1] vmax = rangevals[const][gain][1]
else: else:
vmin = np.nanmin(np.array(mod_data[key])) vmin = np.nanmin(np.array(mod_data[key]))
vmax = np.nanmean( vmax = np.nanmean(
np.array(mod_data[key])) + 2*np.nanstd(np.array(mod_data[key])) np.array(mod_data[key])) + 2*np.nanstd(np.array(mod_data[key]))
if const in ['SlopesFF', 'SlopesPC', 'SlopesCI']: if const in ['SlopesFF', 'SlopesPC', 'SlopesCI']:
htype = 2 htype = 2
else: else:
htype = 0 htype = 0
if key == 'NBP': if key == 'NBP':
unit = '[%]' unit = '[%]'
else: else:
unit = '[ADU]' unit = '[ADU]'
if const == 'Noise-e': if const == 'Noise-e':
unit = '[$e^-$]' unit = '[$e^-$]'
title = '{}, All modules, {} gain, {}'.format( title = '{}, All modules, {} gain, {}'.format(
const, gain_name[gain], keys[key][1]) const, gain_name[gain], keys[key][1])
cb_label = '{}, {} {}'.format(const, keys[key][2], unit) cb_label = '{}, {} {}'.format(const, keys[key][2], unit)
HMCombine(np.array(mod_data[key])[mod_idx][::-1], hm_combine(np.array(mod_data[key])[mod_idx][::-1],
y_ticks=np.arange(nModules)[::-1]+0.8, y_ticks=np.arange(nModules)[::-1]+0.8,
y_ticklabels=np.array(mod_names)[mod_idx], y_ticklabels=np.array(mod_names)[mod_idx],
x_label='Creation Time', y_label='Module ID', x_label='Creation Time', y_label='Module ID',
x_ticklabels=ctimes_ticks, x_ticks=np.arange( x_ticklabels=ctimes_ticks, x_ticks=np.arange(
len(ctimes_ticks))+0.3, len(ctimes_ticks))+0.3,
title=title, cb_label=cb_label, title=title, cb_label=cb_label,
fname='{}/{}_all_g{}_{}.png'.format( fname='{}/{}_all_g{}_{}.png'.format(
out_folder, const, gain, key), out_folder, const, gain, key),
vmin=vmin, vmax=vmax, vmin=vmin, vmax=vmax,
pad=[0.125, 0.151, 0.12, 0.17], type=htype) pad=[0.125, 0.151, 0.12, 0.17], type=htype)
#plt.show()
#break
#break
#break
``` ```
......
notebooks/LPD/PlotFromCalDB_LPD_NBC.ipynb
+ 48
− 60
View file @ bff91e27
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Statistical analysis of calibration factors# # Statistical analysis of calibration factors#
Author: Mikhail Karnevskiy, Steffen Hauf, Version 0.2 Author: Mikhail Karnevskiy, Steffen Hauf, Version 0.2
Calibration constants for LPD1M1 detector from the data base with injection time between start_date and end_date are considered. Calibration constants for LPD1M1 detector from the data base with injection time between start_date and end_date are considered.
To be visualized, calibration constants are averaged per ASICs. Plots shows calibration constant over time for each constant and for each module. Summary plots overall modules are created. To be visualized, calibration constants are averaged per ASICs. Plots shows calibration constant over time for each constant and for each module. Summary plots overall modules are created.
In additional gain-slopes flat-field and pulse-capacitor are combined to relative-gain constant and presented as well. Noise in electron units is derived using gain factors and presented. In additional gain-slopes flat-field and pulse-capacitor are combined to relative-gain constant and presented as well. Noise in electron units is derived using gain factors and presented.
Values shown in plots are saved in h5 files. Values shown in plots are saved in h5 files.
All presented values corresponds to high and medium gain stages. All presented values corresponds to high and medium gain stages.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
cluster_profile = "noDB" # The ipcluster profile to use cluster_profile = "noDB" # The ipcluster profile to use
start_date = "2018-01-30" # Date to start investigation interval from start_date = "2018-01-30" # Date to start investigation interval from
end_date = "2018-12-12" # Date to end investigation interval at, can be "now" end_date = "2018-12-12" # Date to end investigation interval at, can be "now"
constants = ["Offset", "Noise", "SlopesFF", "SlopesCI"] # constants to plot constants = ["Offset", "Noise", "SlopesFF", "SlopesCI"] # constants to plot
modules = [1] # Modules, range allowed modules = [2] # Modules, range allowed
bias_voltages = [250, 500] # Bias voltage bias_voltages = [250, 500] # Bias voltage
mem_cells = [1, 128, 256, 512] # Number of used memory cells. Typically: 4,32,64,128,176. mem_cells = [1, 128, 256, 512] # Number of used memory cells. Typically: 4,32,64,128,176.
photon_energy = 9.2 # Photon energy of the beam photon_energy = 9.2 # Photon energy of the beam
out_folder = "/gpfs/exfel/data/scratch/karnem/testLPD_16/" # Output folder, required out_folder = "/gpfs/exfel/data/scratch/karnem/testLPD_15/" # Output folder, required
use_existing = "" # Input folder use_existing = "/gpfs/exfel/data/scratch/karnem/testLPD_16/" # Input folder
cal_db_timeout = 180000 # timeout on caldb requests", cal_db_timeout = 180000 # timeout on caldb requests",
detectorDB = "LPD1M1" # detector entry in the DB to investigate detectorDB = "LPD1M1" # detector entry in the DB to investigate
dclass = "LPD" # Detector class dclass = "LPD" # Detector class
cal_db_interface = "tcp://max-exfl016:8015#8025" # the database interface to use cal_db_interface = "tcp://max-exfl016:8015#8025" # the database interface to use
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import os
import numpy as np
import matplotlib
##matplotlib.use("agg")
import matplotlib.pyplot as plt
% matplotlib inline
import sys
import datetime import datetime
import dateutil.parser import dateutil.parser
import numpy as np
import os
import sys
import warnings import warnings
warnings.filterwarnings('ignore') warnings.filterwarnings('ignore')
import matplotlib
% matplotlib inline
from iCalibrationDB import Constants, Conditions, Detectors from iCalibrationDB import Constants, Conditions, Detectors
from cal_tools.tools import get_from_db from cal_tools.tools import get_from_db
from cal_tools.ana_tools import * from cal_tools.ana_tools import *
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Prepare variables # Prepare variables
interval = 1 # interval for evaluation in days interval = 1 # interval for evaluation in days
if modules[0] == -1: nMem = 512 # Number of mem Cells to store
modules = list(range(16)) spShape = (64,64) # Shape of superpixel
adu_to_photon = 8000 / 3.6 / 67.# ADU to photon conversion factor
in_mod_names = [] in_mod_names = []
for mod in modules: for mod in modules:
in_mod_names.append("Q{}M{}".format(mod // 4 + 1, mod % 4 + 1)) in_mod_names.append("Q{}M{}".format(mod // 4 + 1, mod % 4 + 1))
constantsDark = {"SlopesFF": 'BadPixelsFF', constantsDark = {"SlopesFF": 'BadPixelsFF',
'SlopesCI': 'BadPixelsCI', 'SlopesCI': 'BadPixelsCI',
'Noise': 'BadPixelsDark', 'Noise': 'BadPixelsDark',
'Offset': 'BadPixelsDark'} 'Offset': 'BadPixelsDark'}
print('Bad pixels data: ', constantsDark) print('Bad pixels data: ', constantsDark)
# Define parameters in order to perform loop over time stamps # Define parameters in order to perform loop over time stamps
start = datetime.datetime.now() if start_date.upper() == "NOW" else dateutil.parser.parse( start = datetime.datetime.now() if start_date.upper() == "NOW" else dateutil.parser.parse(
start_date) start_date)
end = datetime.datetime.now() if end_date.upper() == "NOW" else dateutil.parser.parse( end = datetime.datetime.now() if end_date.upper() == "NOW" else dateutil.parser.parse(
end_date) end_date)
step = datetime.timedelta(hours=interval) step = datetime.timedelta(hours=interval)
dt = end dt = end
# Create output folder # Create output folder
os.makedirs(out_folder, exist_ok=True) os.makedirs(out_folder, exist_ok=True)
# Get getector conditions # Get getector conditions
det = getattr(Detectors, detectorDB) det = getattr(Detectors, detectorDB)
dconstants = getattr(Constants, dclass) dconstants = getattr(Constants, dclass)
if "#" in cal_db_interface: if "#" in cal_db_interface:
prot, serv, ran = cal_db_interface.split(":") prot, serv, ran = cal_db_interface.split(":")
r1, r2 = ran.split("#") r1, r2 = ran.split("#")
cal_db_interface = ":".join( cal_db_interface = ":".join(
[prot, serv, str(np.random.randint(int(r1), int(r2)))]) [prot, serv, str(np.random.randint(int(r1), int(r2)))])
print('CalDB Interface {}'.format(cal_db_interface)) print('CalDB Interface {}'.format(cal_db_interface))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import copy import copy
def prepare_to_store(a, nMem): def prepare_to_store(a, nMem):
shape = list(a.shape[:2])+[nMem,2] shape = list(a.shape[:2])+[nMem, 2]
b = np.full(shape, np.nan) b = np.full(shape, np.nan)
b[:, :, :a.shape[2]] = a[:, :, :, :2] b[:, :, :a.shape[2]] = a[:, :, :, :2]
return b return b
def get_rebined(a, rebin): def get_rebined(a, rebin):
return a.reshape( return a.reshape(
int(a.shape[0] / rebin[0]), int(a.shape[0] / rebin[0]),
rebin[0], rebin[0],
int(a.shape[1] / rebin[1]), int(a.shape[1] / rebin[1]),
rebin[1], rebin[1],
a.shape[2], a.shape[2],
a.shape[3]) a.shape[3])
nMem = 512 # Number of mem Cells to store def modify_const(const, data):
spShape = (64,64) # Shape of superpixel
if const in ['SlopesFF']:
data = data[..., None, None]
if(len(data.shape)==5):
data = data[:,:,:,:,0]
if len(data.shape) < 4:
print(data.shape, "Unexpected shape !!!!!!!!!")
if data.shape[0] != 256:
data = data.swapaxes(0, 2).swapaxes(1,3).swapaxes(2,3)
return data
ret_constants = {} ret_constants = {}
if use_existing != '': if use_existing != '':
mem_cells = [] mem_cells = []
# Loop over max_memory cells # Loop over max_memory cells
for the_mem_cells in mem_cells: for the_mem_cells in mem_cells:
# Loop over bias voltages # Loop over bias voltages
for bias_voltage in bias_voltages: for bias_voltage in bias_voltages:
# Loop over constants # Loop over constants
for c, const in enumerate(constants): for c, const in enumerate(constants):
if not const in ret_constants: if not const in ret_constants:
ret_constants[const] = {} ret_constants[const] = {}
# Loop over modules # Loop over modules
for module in modules: for module in modules:
dt = end dt = end
qm = "Q{}M{}".format(module // 4 + 1, module % 4 + 1) qm = "Q{}M{}".format(module // 4 + 1, module % 4 + 1)
if not qm in ret_constants[const]: if not qm in ret_constants[const]:
ret_constants[const][qm] = [] ret_constants[const][qm] = []
# Loop over time stamps # Loop over time stamps
while dt > start: while dt > start:
creation_time = dt creation_time = dt
if (const in ["Offset", "Noise", if (const in ["Offset", "Noise",
"SlopesCI"] or "DARK" in const.upper()): "SlopesCI"] or "DARK" in const.upper()):
dcond = Conditions.Dark dcond = Conditions.Dark
mcond = getattr(dcond, dclass)( mcond = getattr(dcond, dclass)(
memory_cells=the_mem_cells, memory_cells=the_mem_cells,
bias_voltage=bias_voltage) bias_voltage=bias_voltage)
else: else:
dcond = Conditions.Illuminated dcond = Conditions.Illuminated
mcond = getattr(dcond, dclass)( mcond = getattr(dcond, dclass)(
memory_cells=the_mem_cells, memory_cells=the_mem_cells,
bias_voltage=bias_voltage, bias_voltage=bias_voltage,
photon_energy=photon_energy) photon_energy=photon_energy)
print('Request: ', const, qm, the_mem_cells, bias_voltage, print('Request: ', const, qm, the_mem_cells, bias_voltage,
creation_time) creation_time)
cdata, ctime = get_from_db(getattr(det, qm), cdata, ctime = get_from_db(getattr(det, qm),
getattr(dconstants, const)(), getattr(dconstants, const)(),
copy.deepcopy(mcond), copy.deepcopy(mcond),
None, None,
cal_db_interface, cal_db_interface,
creation_time=creation_time, creation_time=creation_time,
verbosity=0, verbosity=0,
timeout=cal_db_timeout, timeout=cal_db_timeout,
meta_only=True ) meta_only=True )
if ctime is None or cdata is None: if ctime is None or cdata is None:
print('Time or Data is None') print('Time or Data is None')
break break
ctime = ctime.calibration_constant_version.begin_at ctime = ctime.calibration_constant_version.begin_at
print("Found constant {}: begin at {}".format(const, print("Found constant {}: begin at {}".format(const,
cdata is not None), cdata is not None),
ctime) ctime)
print(cdata.shape) print(cdata.shape)
cdata = modify_const(const, cdata)
if const in ['SlopesFF']:
cdata = cdata[...,None,None]
if(len(cdata.shape)==5):
cdata = cdata[:,:,:,:,0]
if len(cdata.shape) < 4:
print(cdata.shape, "Unexpected shape !!!!!!!!!")
if cdata.shape[0] != 256:
cdata = cdata.swapaxes(0, 2).swapaxes(1,3).swapaxes(2,3)
print(cdata.shape) print(cdata.shape)
# Request bad pixel mask # Request bad pixel mask
if const in constantsDark: if const in constantsDark:
#for par in mcond.parameters:
# print (par.name, par.value)
print('RequestBP ', print('RequestBP ',
(creation_time + step + step + step + step)) (creation_time + step + step + step + step))
cdataBP, ctimeBP = get_from_db(getattr(det, qm), cdataBP, ctimeBP = get_from_db(getattr(det, qm),
getattr(dconstants, getattr(dconstants,
constantsDark[const])(), constantsDark[const])(),
copy.deepcopy(mcond), copy.deepcopy(mcond),
None, None,
cal_db_interface, cal_db_interface,
creation_time=( creation_time=(
creation_time + step + step + step + step), creation_time + step + step + step + step),
verbosity=0, verbosity=0,
timeout=cal_db_timeout, timeout=cal_db_timeout,
meta_only=True) meta_only=True)
ctimeBP = ctimeBP.calibration_constant_version.begin_at ctimeBP = ctimeBP.calibration_constant_version.begin_at
if cdataBP is None: if cdataBP is None:
print("Bad pixels are not found !!!!!!!!!") print("Bad pixels are not found !!!!!!!!!")
ctime = ctime.replace(tzinfo=None) ctime = ctime.replace(tzinfo=None)
dt = ctime - step dt = ctime - step
continue continue
print("Found BP {}".format(ctimeBP)) print("Found BP {}".format(ctimeBP))
print(cdataBP.shape) print(cdataBP.shape)
cdataBP = modify_const(const, cdataBP)
if const in ['SlopesFF']: print(cdataBP.shape)
cdataBP = cdataBP[...,None,None]
if(len(cdataBP.shape)==5):
cdataBP = cdataBP[:,:,:,:,0]
if cdataBP.shape[0] != 256:
cdataBP = cdataBP.swapaxes(0, 2).swapaxes(1,3).swapaxes(2,3)
if (len(cdataBP.shape) < 4):
print(cdataBP.shape, "Unexpected shape !!!!!!!!!")
if cdataBP.shape != cdata.shape: if cdataBP.shape != cdata.shape:
print(cdataBP.shape, 'Wrong BP shape!!!') print('Wrong BP shape!!!')
ctime = ctime.replace(tzinfo=None) ctime = ctime.replace(tzinfo=None)
dt = ctime - step dt = ctime - step
continue continue
# Apply bad pixel mask # Apply bad pixel mask
cdataABP = np.copy(cdata) cdataABP = np.copy(cdata)
cdataABP[cdataBP > 0] = np.nan cdataABP[cdataBP > 0] = np.nan
# Create superpixels for constants with BP applied # Create superpixels for constants with BP applied
cdataABP = get_rebined(cdataABP, spShape) cdataABP = get_rebined(cdataABP, spShape)
toStoreBP = prepare_to_store(np.nanmean(cdataABP, axis=(1, 3)), nMem) toStoreBP = prepare_to_store(np.nanmean(cdataABP, axis=(1, 3)), nMem)
toStoreBPStd = prepare_to_store(np.nanstd(cdataABP, axis=(1, 3)), nMem) toStoreBPStd = prepare_to_store(np.nanstd(cdataABP, axis=(1, 3)), nMem)
# Prepare number of bad pixels per superpixels # Prepare number of bad pixels per superpixels
cdataBP = get_rebined(cdataBP, spShape) cdataBP = get_rebined(cdataBP, spShape)
cdataNBP = prepare_to_store(np.nansum(cdataBP > 0, axis=(1, 3)), nMem) cdataNBP = prepare_to_store(np.nansum(cdataBP > 0, axis=(1, 3)), nMem)
else: else:
toStoreExtBP = 0 toStoreExtBP = 0
cdataBPExt = 0 cdataBPExt = 0
# Create superpixels for constants without BP applied # Create superpixels for constants without BP applied
cdata = get_rebined(cdata, spShape) cdata = get_rebined(cdata, spShape)
toStoreStd = prepare_to_store(np.nanstd(cdata, axis=(1, 3)), nMem) toStoreStd = prepare_to_store(np.nanstd(cdata, axis=(1, 3)), nMem)
toStore = prepare_to_store(np.nanmean(cdata, axis=(1, 3)), nMem) toStore = prepare_to_store(np.nanmean(cdata, axis=(1, 3)), nMem)
# Convert parameters to dict # Convert parameters to dict
dpar = {} dpar = {}
for par in mcond.parameters: for par in mcond.parameters:
dpar[par.name] = par.value dpar[par.name] = par.value
print("Store values in dict", const, qm, ctime) print("Store values in dict", const, qm, ctime)
ret_constants[const][qm].append({'ctime': ctime, ret_constants[const][qm].append({'ctime': ctime,
'nBP': cdataNBP, 'nBP': cdataNBP,
'dataBP': toStoreBP, 'dataBP': toStoreBP,
'dataBPStd': toStoreBPStd, 'dataBPStd': toStoreBPStd,
'data': toStore, 'data': toStore,
'dataStd': toStoreStd, 'dataStd': toStoreStd,
'mdata': dpar}) 'mdata': dpar})
ctime = ctime.replace(tzinfo=None) ctime = ctime.replace(tzinfo=None)
dt = ctime - step dt = ctime - step
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if use_existing == "": if use_existing == "":
print('Save data to /CalDBAna_{}_{}.h5'.format(dclass, modules[0])) print('Save data to /CalDBAna_{}_{}.h5'.format(dclass, modules[0]))
save_dict_to_hdf5(ret_constants, save_dict_to_hdf5(ret_constants,
'{}/CalDBAna_{}_{}.h5'.format(out_folder, dclass, modules[0])) '{}/CalDBAna_{}_{}.h5'.format(out_folder, dclass, modules[0]))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if use_existing == "": if use_existing == "":
fpath = '{}/CalDBAna_{}_*.h5'.format(out_folder, dclass) fpath = '{}/CalDBAna_{}_*.h5'.format(out_folder, dclass)
else: else:
fpath = '{}/CalDBAna_{}_*.h5'.format(use_existing, dclass) fpath = '{}/CalDBAna_{}_*.h5'.format(use_existing, dclass)
print('Load data from {}'.format(fpath)) print('Load data from {}'.format(fpath))
ret_constants = load_data_from_hdf5(fpath) ret_constants = load_data_from_hdf5(fpath)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
print ('Combine Gain FF and CI and Noise in [e-]') # Combine FF and PC data to calculate Gain
# Estimate Noise in units of electrons
print ('Calculate Gain and Noise in electron units')
ret_constants["Gain"] = {} ret_constants["Gain"] = {}
ret_constants["Noise-e"] = {} ret_constants["Noise-e"] = {}
for module in list(range(16)): for module in list(range(16)):
if ("SlopesFF" not in ret_constants or if ("SlopesFF" not in ret_constants or
"SlopesPC" not in ret_constants): "SlopesCI" not in ret_constants):
break break
qm = "Q{}M{}".format(module // 4 + 1, module % 4 + 1) qm = "Q{}M{}".format(module // 4 + 1, module % 4 + 1)
print(qm) print(qm)
if (qm not in ret_constants["SlopesFF"] or if (qm not in ret_constants["SlopesFF"] or
qm not in ret_constants["SlopesCI"]): qm not in ret_constants["SlopesCI"]):
continue continue
ret_constants["Gain"][qm] = {} ret_constants["Gain"][qm] = {}
dataFF = ret_constants["SlopesFF"][qm] dataFF = ret_constants["SlopesFF"][qm]
dataPC = ret_constants["SlopesCI"][qm] dataPC = ret_constants["SlopesCI"][qm]
if (len(dataFF) == 0 or len(dataPC) == 0): if (len(dataFF) == 0 or len(dataPC) == 0):
continue continue
ctimesFF = np.array(dataFF["ctime"]) ctimesFF = np.array(dataFF["ctime"])
ctimesPC = np.array(dataPC["ctime"]) ctimesPC = np.array(dataPC["ctime"])
ctime, icomb = combine_constants(ctimesFF, ctimesPC) ctime, icomb = combine_constants(ctimesFF, ctimesPC)
cdataPC_vs_time = np.array(dataPC["data"])[..., 0] cdataPC_vs_time = np.array(dataPC["data"])[..., 0]
cdataFF_vs_time = np.array(dataFF["data"])[..., 0] cdataFF_vs_time = np.array(dataFF["data"])[..., 0]
cdataFF_vs_time = np.nanmedian(cdataFF_vs_time, axis=3)[..., None] cdataFF_vs_time = np.nanmedian(cdataFF_vs_time, axis=3)[..., None]
cdataFF_vs_time /= np.nanmedian(cdataFF_vs_time, axis=(1, 2, 3))[:, None, cdataFF_vs_time /= np.nanmedian(cdataFF_vs_time, axis=(1, 2, 3))[:, None,
None, None] None, None]
cdataPC_vs_time /= np.nanmedian(cdataPC_vs_time, axis=(1, 2, 3))[:, None, cdataPC_vs_time /= np.nanmedian(cdataPC_vs_time, axis=(1, 2, 3))[:, None,
None, None] None, None]
gain_vs_time = [] gain_vs_time = []
for iFF, iPC in icomb: for iFF, iPC in icomb:
gain_vs_time.append(cdataFF_vs_time[iFF] * cdataPC_vs_time[iPC]) gain_vs_time.append(cdataFF_vs_time[iFF] * cdataPC_vs_time[iPC])
print(np.array(gain_vs_time).shape) print(np.array(gain_vs_time).shape)
ctime_ts = [t.timestamp() for t in ctime] ctime_ts = [t.timestamp() for t in ctime]
ret_constants["Gain"][qm]["ctime"] = ctime ret_constants["Gain"][qm]["ctime"] = ctime
ret_constants["Gain"][qm]["data"] = np.array(gain_vs_time) ret_constants["Gain"][qm]["data"] = np.array(gain_vs_time)
# Fill missing data for compatibility with plotting code
ret_constants["Gain"][qm]["dataBP"] = np.array(gain_vs_time) ret_constants["Gain"][qm]["dataBP"] = np.array(gain_vs_time)
ret_constants["Gain"][qm]["nBP"] = np.array(gain_vs_time) ret_constants["Gain"][qm]["nBP"] = np.array(gain_vs_time)
if "Noise" not in ret_constants: if "Noise" not in ret_constants:
continue continue
if qm not in ret_constants["Noise"]: if qm not in ret_constants["Noise"]:
continue continue
dataN = ret_constants["Noise"][qm] dataN = ret_constants["Noise"][qm]
if len(dataN) == 0: if len(dataN) == 0:
continue continue
ret_constants["Noise-e"][qm] = {} ret_constants["Noise-e"][qm] = {}
ctimesG = np.array(ctime) ctimesG = np.array(ctime)
ctimesN = np.array(dataN["ctime"]) ctimesN = np.array(dataN["ctime"])
ctime, icomb = combine_constants(ctimesG, ctimesN) ctime, icomb = combine_constants(ctimesG, ctimesN)
cdataG_vs_time = np.array(gain_vs_time) cdataG_vs_time = np.array(gain_vs_time)
cdataN_vs_time = np.array(dataN["data"])[..., 0] cdataN_vs_time = np.array(dataN["data"])[..., 0]
data_vs_time = [] data_vs_time = []
for iG, iN in icomb: for iG, iN in icomb:
data_vs_time.append( data_vs_time.append(
cdataN_vs_time[iN] * (8000 / 3.6 / 67.) / cdataG_vs_time[iG]) cdataN_vs_time[iN] * adu_to_photon / cdataG_vs_time[iG])
print(np.array(gain_vs_time).shape) print(np.array(gain_vs_time).shape)
ctime_ts = [t.timestamp() for t in ctime] ctime_ts = [t.timestamp() for t in ctime]
ret_constants["Noise-e"][qm]["ctime"] = ctime ret_constants["Noise-e"][qm]["ctime"] = ctime
ret_constants["Noise-e"][qm]["data"] = np.array(data_vs_time) ret_constants["Noise-e"][qm]["data"] = np.array(data_vs_time)
# Fill missing data for compatibility with plotting code
ret_constants["Noise-e"][qm]["dataBP"] = np.array(data_vs_time) ret_constants["Noise-e"][qm]["dataBP"] = np.array(data_vs_time)
ret_constants["Noise-e"][qm]["nBP"] = np.array(data_vs_time) ret_constants["Noise-e"][qm]["nBP"] = np.array(data_vs_time)
save_dict_to_hdf5({k:v for k,v in ret_constants.items() if k in ['Gain', 'Noise-e']}, save_dict_to_hdf5({k:v for k,v in ret_constants.items() if k in ['Gain', 'Noise-e']},
'{}/CalDBAna_{}_Gain.h5'.format(out_folder, dclass)) '{}/CalDBAna_{}_Gain.h5'.format(out_folder, dclass))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Parameters for plotting # Parameters for plotting
# Define range for plotting # Define range for plotting
rangevals = { rangevals = {
"Offset": [[800., 1500], [600, 900]], "Offset": [[800., 1500], [600, 900]],
"Noise": [[2.0, 16], [1.0, 7.0]], "Noise": [[2.0, 16], [1.0, 7.0]],
"Gain": [[20, 30], [20, 30]], "Gain": [[20, 30], [20, 30]],
"Noise-e": [[100., 600.], [100., 600.]], "Noise-e": [[100., 600.], [100., 600.]],
"SlopesCI": [[0.95, 1.05], [0.0, 0.5]], "SlopesCI": [[0.95, 1.05], [0.0, 0.5]],
"SlopesFF": [[0.8, 1.2], [0.8, 1.2]] "SlopesFF": [[0.8, 1.2], [0.8, 1.2]]
} }
nMemToShow = 32 nMemToShow = 32
keys = { keys = {
'Mean': ['data', '', 'Mean over pixels'], 'Mean': ['data', '', 'Mean over pixels'],
'std': ['dataStd', '', '$\sigma$ over pixels'], 'std': ['dataStd', '', '$\sigma$ over pixels'],
'MeanBP': ['dataBP', 'Good pixels only', 'Mean over pixels'], 'MeanBP': ['dataBP', 'Good pixels only', 'Mean over pixels'],
'NBP': ['nBP', 'Fraction of BP', 'Fraction of BP'], 'NBP': ['nBP', 'Fraction of BP', 'Fraction of BP'],
'stdBP': ['dataBPStd', 'Good pixels only', '$\sigma$ over pixels'], 'stdBP': ['dataBPStd', 'Good pixels only', '$\sigma$ over pixels'],
'stdASIC': ['', '', '$\sigma$ over ASICs'], 'stdASIC': ['', '', '$\sigma$ over ASICs'],
'stdCell': ['', '', '$\sigma$ over Cells'], 'stdCell': ['', '', '$\sigma$ over Cells'],
} }
gain_name = ['High', 'Medium', 'Low'] gain_name = ['High', 'Medium', 'Low']
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
print('Plot calibration constants') print('Plot calibration constants')
# loop over constat type # loop over constat type
for const, modules in ret_constants.items(): for const, modules in ret_constants.items():
# Loop over gain # Loop over gain
for gain in range(2): for gain in range(2):
print('Const: {}, gain {}'.format(const, gain)) print('Const: {}, gain {}'.format(const, gain))
if const in ["Gain", "Noise-e"] and gain == 1: if const in ["Gain", "Noise-e"] and gain == 1:
continue continue
else: else:
pass pass
# loop over modules # loop over modules
mod_data = {} mod_data = {}
mod_data['stdASIC'] = [] mod_data['stdASIC'] = []
mod_data['stdCell'] = [] mod_data['stdCell'] = []
mod_names = [] mod_names = []
mod_times = [] mod_times = []
# Loop over modules # Loop over modules
for mod, data in modules.items(): for mod, data in modules.items():
if mod not in in_mod_names: if mod not in in_mod_names:
continue continue
print(mod) print(mod)
ctimes = np.array(data["ctime"]) ctimes = np.array(data["ctime"])
ctimes_ticks = [x.strftime('%m-%d') for x in ctimes] ctimes_ticks = [x.strftime('%m-%d') for x in ctimes]
if ("mdata" in data): if ("mdata" in data):
cmdata = np.array(data["mdata"]) cmdata = np.array(data["mdata"])
for i, tick in enumerate(ctimes_ticks): for i, tick in enumerate(ctimes_ticks):
ctimes_ticks[i] = ctimes_ticks[i] + \ ctimes_ticks[i] = ctimes_ticks[i] + \
', V={:1.0f}'.format(cmdata[i]['Sensor Bias Voltage']) + \ ', V={:1.0f}'.format(cmdata[i]['Sensor Bias Voltage']) + \
', M={:1.0f}'.format( ', M={:1.0f}'.format(
cmdata[i]['Memory cells']) cmdata[i]['Memory cells'])
sort_ind = np.argsort(ctimes_ticks) sort_ind = np.argsort(ctimes_ticks)
ctimes_ticks = list(np.array(ctimes_ticks)[sort_ind]) ctimes_ticks = list(np.array(ctimes_ticks)[sort_ind])
# Create sorted by data dataset # Create sorted by data dataset
rdata = {} rdata = {}
for key, item in keys.items(): for key, item in keys.items():
if item[0] in data: if item[0] in data:
rdata[key] = np.array(data[item[0]])[sort_ind] rdata[key] = np.array(data[item[0]])[sort_ind]
# print(rdata['Mean'].shape)
nTimes = rdata['Mean'].shape[0] nTimes = rdata['Mean'].shape[0]
nPixels = rdata['Mean'].shape[1] * rdata['Mean'].shape[2] nPixels = rdata['Mean'].shape[1] * rdata['Mean'].shape[2]
nBins = nMemToShow * nPixels nBins = nMemToShow * nPixels
# Select gain # Select gain
if const not in ["Gain", "Noise-e"]: if const not in ["Gain", "Noise-e"]:
for key in rdata: for key in rdata:
rdata[key] = rdata[key][..., gain] rdata[key] = rdata[key][..., gain]
# Avoid to low values # Avoid to low values
if const in ["Noise", "Offset", "Noise-e"]: if const in ["Noise", "Offset", "Noise-e"]:
rdata['Mean'][rdata['Mean'] < 0.1] = np.nan rdata['Mean'][rdata['Mean'] < 0.1] = np.nan
if 'MeanBP' in rdata: if 'MeanBP' in rdata:
rdata['MeanBP'][rdata['MeanBP'] < 0.1] = np.nan rdata['MeanBP'][rdata['MeanBP'] < 0.1] = np.nan
if 'NBP' in rdata: if 'NBP' in rdata:
rdata["NBP"][rdata["NBP"] == 4096] = np.nan rdata["NBP"][rdata["NBP"] == 4096] = np.nan
rdata["NBP"] = rdata["NBP"] / (64 * 64) * 100 rdata["NBP"] = rdata["NBP"] / (64 * 64) * 100
# Reshape: ASICs over cells for plotting # Reshape: ASICs over cells for plotting
pdata = {} pdata = {}
for key in rdata: for key in rdata:
pdata[key] = rdata[key][:, :, :, :nMemToShow].reshape( pdata[key] = rdata[key][:, :, :, :nMemToShow].reshape(
nTimes, nBins).swapaxes(0, 1) nTimes, nBins).swapaxes(0, 1)
# Summary over ASICs # Summary over ASICs
adata = {} adata = {}
for key in rdata: for key in rdata:
adata[key] = np.nanmean(rdata[key], axis=(1, 2)).swapaxes(0, 1) adata[key] = np.nanmean(rdata[key], axis=(1, 2)).swapaxes(0, 1)
# Summary information over modules # Summary information over modules
for key in pdata: for key in pdata:
if key not in mod_data: if key not in mod_data:
mod_data[key] = [] mod_data[key] = []
mod_data[key].append(np.nanmean(pdata[key], axis=0)) mod_data[key].append(np.nanmean(pdata[key], axis=0))
mod_data['stdASIC'].append(np.nanstd( mod_data['stdASIC'].append(np.nanstd(
np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=(1, 2)), axis=1)) np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=(1, 2)), axis=1))
mod_data['stdCell'].append(np.nanstd( mod_data['stdCell'].append(np.nanstd(
np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=3), axis=(1, 2))) np.nanmean(rdata['Mean'][:, :, :, :nMemToShow], axis=3), axis=(1, 2)))
mod_names.append(mod) mod_names.append(mod)
mod_times.append(ctimes_ticks) mod_times.append(ctimes_ticks)
# Plotting # Plotting
for key in pdata: for key in pdata:
vmin = None vmin = None
vmax = None vmax = None
if const in rangevals and key in ['Mean', 'MeanBP']: if const in rangevals and key in ['Mean', 'MeanBP']:
vmin = rangevals[const][gain][0] vmin = rangevals[const][gain][0]
vmax = rangevals[const][gain][1] vmax = rangevals[const][gain][1]
if key == 'NBP': if key == 'NBP':
unit = '[%]' unit = '[%]'
else: else:
unit = '[ADU]' unit = '[ADU]'
if const == 'Noise-e': if const == 'Noise-e':
unit = '[$e^-$]' unit = '[$e^-$]'
title = '{}, module {}, {} gain, {}'.format( title = '{}, module {}, {} gain, {}'.format(
const, mod, gain_name[gain], keys[key][1]) const, mod, gain_name[gain], keys[key][1])
cb_label = '{}, {} {}'.format(const, keys[key][2], unit) cb_label = '{}, {} {}'.format(const, keys[key][2], unit)
HMCombine(pdata[key][::-1], HMCombine(pdata[key][::-1],
x_label='Creation Time', y_label='ASIC ID', x_label='Creation Time', y_label='ASIC ID',
x_ticklabels=ctimes_ticks, x_ticklabels=ctimes_ticks,
x_ticks=np.arange(len(ctimes_ticks))+0.3, x_ticks=np.arange(len(ctimes_ticks))+0.3,
title=title, cb_label=cb_label, title=title, cb_label=cb_label,
vmin=vmin, vmax=vmax, vmin=vmin, vmax=vmax,
fname='{}/{}_{}_g{}_ASIC_{}.png'.format( fname='{}/{}_{}_g{}_ASIC_{}.png'.format(
out_folder, const, mod, gain, key), out_folder, const, mod, gain, key),
y_ticks=np.arange(nBins, step=nMemToShow)+16, y_ticks=np.arange(nBins, step=nMemToShow)+16,
y_ticklabels=np.arange(nPixels)[::-1]+1, y_ticklabels=np.arange(nPixels)[::-1]+1,
pad=[0.125, 0.125, 0.12, 0.185]) pad=[0.125, 0.125, 0.12, 0.185])
HMCombine(adata[key], type=0, HMCombine(adata[key], type=0,
x_label='Creation Time', y_label='Memory cell ID', x_label='Creation Time', y_label='Memory cell ID',
x_ticklabels=ctimes_ticks, x_ticklabels=ctimes_ticks,
x_ticks=np.arange(len(ctimes_ticks))+0.3, x_ticks=np.arange(len(ctimes_ticks))+0.3,
title=title, cb_label=cb_label, title=title, cb_label=cb_label,
fname='{}/{}_{}_g{}_MEM_{}.png'.format( fname='{}/{}_{}_g{}_MEM_{}.png'.format(
out_folder, const, mod, gain, key), out_folder, const, mod, gain, key),
vmin=vmin, vmax=vmax) vmin=vmin, vmax=vmax)
#break
#break
#break
#break
``` ```
......
xfel_calibrate/notebooks.py
+ 10
− 2
View file @ bff91e27
...@@ -33,6 +33,13 @@ notebooks = { ...@@ -33,6 +33,13 @@ notebooks = {
"default concurrency": None, "default concurrency": None,
"cluster cores": 8}, "cluster cores": 8},
}, },
"STATS_FROM_DB": {
"notebook": "notebooks/AGIPD/PlotFromCalDB_AGIPD_NBC.ipynb",
"dep_notebooks": ["notebooks/AGIPD/PlotFromCalDB_Summary_AGIPD_NBC.ipynb"],
"concurrency": {"parameter": "modules",
"default concurrency": None,
"cluster cores": 1},
},
}, },
"LPD": { "LPD": {
"DARK": { "DARK": {
...@@ -68,8 +75,9 @@ notebooks = { ...@@ -68,8 +75,9 @@ notebooks = {
"cluster cores": 2}, "cluster cores": 2},
}, },
"STATS_FROM_DB": { "STATS_FROM_DB": {
"notebook": "notebooks/LPD/PlotFromCalDB.ipynb", "notebook": "notebooks/LPD/PlotFromCalDB_LPD_NBC.ipynb",
"concurrency": {"parameter": None, "dep_notebooks": ["notebooks/AGIPD/PlotFromCalDB_Summary_AGIPD_NBC.ipynb"],
"concurrency": {"parameter": "modules",
"default concurrency": None, "default concurrency": None,
"cluster cores": 1}, "cluster cores": 1},
}, },
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Legal Notice | Contacts