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

Merge branch 'feat/jungfr_use_karabo_da' into 'master' 

Feat: Jungfrau update: use karabo-da and karabo_id

See merge request detectors/pycalibration!279
parents f578f657 f01ae0c8
No related branches found
No related tags found
2 merge requests!298Feat/dss cimprove master rebasing,!279Feat: Jungfrau update: use karabo-da and karabo_id
Hide whitespace changes
Inline Side-by-side
notebooks/Jungfrau/Jungfrau_Gain_Correct_and_Verify_NBC.ipynb
+ 39
− 26
View file @ 655e14b1
%% Cell type:markdown id: tags:
# Jungfrau Offline Correction #
Author: European XFEL Detector Group, Version: 0.1
Offline Calibration for the Jungfrau Detector
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/SPB/201922/p002566/raw" # the folder to read data from, required
run = 70 # runs to process, required
cluster_profile = "noDB" # cluster profile to use
in_folder = "/gpfs/exfel/exp/CALLAB/202031/p900113/raw" # the folder to read data from, required
out_folder = "/gpfs/exfel/data/scratch/karnem/test/005" # the folder to output to, required
sequences = [-1] # sequences to correct, set to -1 for all, range allowed
run = 9999 # runs to process, required
karabo_id = "FXE_XAD_JF1M" # karabo prefix of Jungfrau devices
karabo_da = 'JNGFR01' # data aggregators
receiver_id = "RECEIVER-{}" # inset for receiver devices
receiver_control_id = "CONTROL" # inset for control devices
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # template to use for file name, double escape sequence number
h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data' # path in H5 file under which images are located
h5path_run = '/RUN/{}/DET/{}' # path to run data
h5path_cntrl = '/CONTROL/{}/DET/{}' # path to control data
karabo_id_control = "" # if control is on a different ID, set to empty string if it is the same a karabo-id
karabo_da_control = "JNGFR01" # file inset for control data
use_dir_creation_date = True # use the creation data of the input dir for database queries
cal_db_interface = "tcp://max-exfl017:8016" #"tcp://max-exfl016:8015#8025" # the database interface to use
cal_db_timeout = 180000 # timeout on caldb requests",
overwrite = True # set to True if existing data should be overwritten
no_relative_gain = False # do not do relative gain correction
cluster_profile = "noDB" # cluster profile to use
bias_voltage = 180 # will be overwritten by value in file
cal_db_interface = "tcp://max-exfl016:8016" #"tcp://max-exfl016:8015#8025" # the database interface to use
cal_db_timeout = 180000 # timeout on caldb requests",
use_dir_creation_date = True # use the creation data of the input dir for database queries
sequences_per_node = 5 # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel
photon_energy = 9.2 # photon energy in keV
chunk_size_idim = 1 # chunking size of imaging dimension, adjust if user software is sensitive to this.
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # template to use for file name, double escape sequence number
path_template_control = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # template to use for file name, double escape sequence number
integration_time = 4.96 # integration time in us, will be overwritten by value in file
h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data' # path in H5 file under which images are located
gmapfile = "" # variable is not used but left here for back compatibility
karabo_id = "SPB_IRDA_JNGFR" # karabo prefix of Jungfrau devices
karabo_id_control = "" # if control is on a different ID, set to empty string for using the same as for data
receiver_id = "MODULE_3" # inset for receiver devices
control_id = "CONTROL" # inset for control devices
db_module = "Jungfrau_M273" # ID of module in calibration database
path_inset = "JNGFR03" # file inset for image data
path_inset_control = "JNGFR01" # file inset for control data
db_module = ["Jungfrau_M233"] # ID of module in calibration database
manual_slow_data = False # if true, use manually entered bias_voltage and integration_time values
chunk_size = 0
def balance_sequences(in_folder, run, sequences, sequences_per_node, path_inset):
def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da):
from xfel_calibrate.calibrate import balance_sequences as bs
return bs(in_folder, run, sequences, sequences_per_node, path_inset)
return bs(in_folder, run, sequences, sequences_per_node, karabo_da)
```
%% Cell type:code id: tags:
``` python
import warnings
warnings.filterwarnings('ignore')
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import XFELDetAna.xfelpyanatools as xana
import XFELDetAna.xfelpycaltools as xcal
from cal_tools.enums import BadPixels
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from cal_tools.tools import (get_dir_creation_date, get_constant_from_db_and_time)
from matplotlib.colors import LogNorm
import os
import h5py
from h5py import File as h5file
import copy
from IPython.display import HTML, display, Markdown, Latex
import tabulate
from XFELDetAna.util import env
env.iprofile = cluster_profile
from functools import partial
from ipyparallel import Client
```
%% Cell type:code id: tags:
``` python
client = Client(profile=cluster_profile)
view = client[:]
view.use_dill()
receiver_id = receiver_id.format(int(karabo_da[-2:]))
h5path = h5path.format(karabo_id, receiver_id)
ped_dir = "{}/r{:04d}".format(in_folder, run)
# TODO
# this trick is needed until proper mapping is introduced
if len(db_module)>1:
db_module = db_module[int(karabo_da[-2:])-1]
else:
db_module = db_module[0]
if ped_dir[-1] == "/":
ped_dir = ped_dir[:-1]
if not os.path.exists(out_folder):
os.makedirs(out_folder)
elif not overwrite:
raise AttributeError("Output path exists! Exiting")
fp_name = path_template.format(run, path_inset)
fp_name = path_template.format(run, karabo_da)
fp_path = '{}/{}'.format(ped_dir, fp_name)
fp_name_contr = path_template_control.format(run, path_inset_control)
fp_name_contr = path_template.format(run, karabo_da_control)
fp_path_contr = '{}/{}'.format(ped_dir, fp_name_contr)
if sequences[0] == -1:
sequences = None
print("Reading data from {}".format(fp_path))
print("Run is: {}".format(run))
print("HDF5 path: {}".format(h5path))
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run)
print("Using {} as creation time".format(creation_time))
if karabo_id_control == "":
karabo_id_control = karabo_id
```
%% Cell type:code id: tags:
``` python
def check_memoryCells(file_name, path):
with h5file(file_name, 'r') as f:
t_stamp = np.array(f[path + '/storageCells/timestamp'])
st_cells = np.array(f[path + '/storageCells/value'])
sc_start = np.array(f[path + '/storageCellStart/value'])
valid_train = t_stamp > 0
n_scs = st_cells[valid_train][0] + 1
sc_s = sc_start[valid_train][0]
return n_scs, sc_s
```
%% Cell type:code id: tags:
``` python
if not manual_slow_data:
with h5py.File(fp_path_contr.format(0), 'r') as f:
integration_time = float(f['/RUN/{}/DET/{}/exposureTime/value'.format(karabo_id_control, control_id)][()]*1e6)
bias_voltage = int(np.squeeze(f['/RUN/{}/DET/{}/vHighVoltage/value'.format(karabo_id_control, control_id)])[0])
run_path = h5path_run.format(karabo_id_control, receiver_control_id)
integration_time = float(f[f'{run_path}/exposureTime/value'][()]*1e6)
bias_voltage = int(np.squeeze(f[f'{run_path}/vHighVoltage/value'])[0])
control_path = '/CONTROL/{}/DET/{}'.format(karabo_id_control, control_id)
control_path = h5path_cntrl.format(karabo_id_control, receiver_control_id)
this_run_mcells, sc_start = check_memoryCells(fp_path_contr.format(0), control_path)
if this_run_mcells == 1:
memoryCells = 1
print('Dark runs in single cell mode\n storage cell start: {:02d}'.format(sc_start))
else:
memoryCells = 16
print('Dark runs in burst mode\n storage cell start: {:02d}'.format(sc_start))
print("Integration time is {} us".format(integration_time))
print("Bias voltage is {} V".format(bias_voltage))
print("Number of memory cells is {}".format(memoryCells))
```
%% Cell type:code id: tags:
``` python
dirlist = sorted(os.listdir(ped_dir))
file_list = []
total_sequences = 0
fsequences = []
for entry in dirlist:
#only h5 file
abs_entry = "{}/{}".format(ped_dir, entry)
if os.path.isfile(abs_entry) and os.path.splitext(abs_entry)[1] == ".h5":
if sequences is None:
for seq in range(len(dirlist)):
if path_template.format(run, path_inset).format(seq) in abs_entry:
if path_template.format(run, karabo_da).format(seq) in abs_entry:
file_list.append(abs_entry)
total_sequences += 1
fsequences.append(seq)
else:
for seq in sequences:
if path_template.format(run, path_inset).format(seq) in abs_entry:
if path_template.format(run, karabo_da).format(seq) in abs_entry:
file_list.append(os.path.abspath(abs_entry))
total_sequences += 1
fsequences.append(seq)
sequences = fsequences
```
%% Cell type:code id: tags:
``` python
print("Processing a total of {} sequence files".format(total_sequences))
table = []
for k, f in enumerate(file_list):
table.append((k, f))
if len(table) > 0:
md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=["#", "file"])))
```
%% Cell type:code id: tags:
``` python
condition = Conditions.Dark.jungfrau(memory_cells=memoryCells,
bias_voltage=bias_voltage,
integration_time=integration_time)
## offset
offset_map, _ = \
get_constant_from_db_and_time(getattr(Detectors, db_module),
Constants.jungfrau.Offset(),
condition,
np.zeros((1024, 512, 1, 3)),
cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout)
mask, _ = \
get_constant_from_db_and_time(getattr(Detectors, db_module),
Constants.jungfrau.BadPixelsDark(),
condition,
np.zeros((1024, 512, 1, 3)),
cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout)
gain_map, _ = \
get_constant_from_db_and_time(getattr(Detectors, db_module),
Constants.jungfrau.RelativeGain(),
condition,
None,
cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout)
# move from x,y,cell,gain to cell,x,y,gain
offset_map = np.squeeze(offset_map)
mask = np.squeeze(mask)
if memoryCells>1:
offset_map = np.moveaxis(np.moveaxis(offset_map, 0, 2), 0, 2)
mask = np.moveaxis(np.moveaxis(mask, 0, 2), 0, 2)
if gain_map is None:
print("No gain map found")
no_relative_gain = True
else:
if memoryCells>1:
gain_map = np.moveaxis(np.moveaxis(gain_map, 0, 2), 0, 1)
else:
gain_map = np.squeeze(gain_map)
gain_map = np.moveaxis(gain_map, 1, 0)
```
%% Cell type:code id: tags:
``` python
def copy_and_sanitize_non_cal_data(infile, outfile, h5base):
""" Copy and sanitize data in `infile` that is not touched by `correctLPD`
"""
if h5base.startswith("/"):
h5base = h5base[1:]
dont_copy = ["adc",]
dont_copy = [h5base+"/{}".format(do)
for do in dont_copy]
def visitor(k, item):
if k not in dont_copy:
if isinstance(item, h5py.Group):
outfile.create_group(k)
elif isinstance(item, h5py.Dataset):
group = str(k).split("/")
group = "/".join(group[:-1])
infile.copy(k, outfile[group])
infile.visititems(visitor)
```
%% Cell type:code id: tags:
``` python
# Loop over sequences
for k, f in enumerate(file_list):
with h5py.File(f, 'r') as infile:
out_file = "{}/{}".format(out_folder, f.split("/")[-1])
out_file = out_file.replace("RAW", "CORR")
print('Process file: ', f)
with h5py.File(out_file, "w") as ofile:
copy_and_sanitize_non_cal_data(infile, ofile, h5path)
oshape = infile[h5path+"/adc"].shape
print('Data shape: ', oshape)
ddset = ofile.create_dataset(h5path+"/adc",
oshape,
chunks=(chunk_size_idim, 1, oshape[1], oshape[2]),
dtype=np.float32)
mskset = ofile.create_dataset(h5path+"/mask",
oshape,
chunks=(chunk_size_idim, 1, oshape[1], oshape[2]),
dtype=np.uint32,
compression="gzip", compression_opts=1, shuffle=True)
# Correct a chunk of images for offset and gain
def correct_chunk(offset_map, mask, gain_map, memoryCells, no_relative_gain, inp):
import numpy as np
import copy
import h5py
fim_data = None
gim_data = None
rim_data = None
msk_data = None
err = ''
try:
d, g, m, ind, copy_sample = inp
g[g==3] = 2
if copy_sample and ind==0:
if memoryCells==1:
rim_data = np.squeeze(copy.copy(d))
else:
rim_data = np.squeeze(copy.copy(d[:,0,...]))
# Select memory cells
if memoryCells>1:
m[m>16] = 0
offset_map_cell = offset_map[m,...]
mask_cell = mask[m,...]
else:
offset_map_cell = offset_map
mask_cell = mask
# Offset correction
offset = np.choose(g, (offset_map_cell[...,0], offset_map_cell[...,1], offset_map_cell[...,2]))
d -= offset
# Gain correction
if not no_relative_gain:
if memoryCells>1:
gain_map_cell = gain_map[m,...]
else:
gain_map_cell = gain_map
cal = np.choose(g, (gain_map_cell[..., 0], gain_map_cell[..., 1], gain_map_cell[..., 2]))
d /= cal
msk = np.choose(g, (mask_cell[...,0], mask_cell[...,1], mask_cell[...,2]))
# Store sample of data for plotting
if copy_sample and ind==0:
if memoryCells==1:
fim_data = np.squeeze(copy.copy(d))
gim_data = np.squeeze(copy.copy(g))
msk_data = np.squeeze(copy.copy(msk))
else:
fim_data = np.squeeze(copy.copy(d[:,1,...]))
gim_data = np.squeeze(copy.copy(g[:,1,...]))
msk_data = np.squeeze(copy.copy(msk[:,1,...]))
print('All done')
except Exception as e:
print(e)
err = e
return ind, d, msk, rim_data, fim_data, gim_data, msk_data, err
# Run ip Cluster parallelization over chunks of images
inp = []
max_ind = oshape[0]
ind = 0
# If chunk size is not given maximum 12+1 chunks is expected
if chunk_size == 0:
chunk_size = max_ind//12
print('Chunk size: {}'.format(chunk_size))
while ind<max_ind:
d = infile[h5path+"/adc"][ind:ind+chunk_size,...].astype(np.float32)
g = infile[h5path+"/gain"][ind:ind+chunk_size,...]
m = infile[h5path+"/memoryCell"][ind:ind+chunk_size,...]
print('To process: ', d.shape)
inp.append((d,g,m, ind, k==0))
ind += chunk_size
print('Run {} processes'.format(len(inp)) )
p = partial(correct_chunk, offset_map, mask, gain_map, memoryCells, no_relative_gain)
r = view.map_sync(p, inp)
#r = list(map(p, inp))
if k==0:
_,_,_, rim_data, fim_data, gim_data, msk_data, _ = r[0]
for rr in r:
ind, cdata, cmask, _,_,_,_, err = rr
data_size = cdata.shape[0]
ddset[ind:ind+data_size,...] = cdata
mskset[ind:ind+data_size,...] = cmask
if err != '':
print('Error: '.format(err))
```
%% Cell type:code id: tags:
``` python
def do_2d_plot(data, edges, y_axis, x_axis):
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
extent = [np.min(edges[1]), np.max(edges[1]),np.min(edges[0]), np.max(edges[0])]
im = ax.imshow(data[::-1,:], extent=extent, aspect="auto", norm=LogNorm(vmin=1, vmax=np.max(data)))
ax.set_xlabel(x_axis)
ax.set_ylabel(y_axis)
cb = fig.colorbar(im)
cb.set_label("Counts")
```
%% Cell type:code id: tags:
``` python
h, ex, ey = np.histogram2d(rim_data.flatten(), gim_data.flatten(),
bins=[100, 4], range=[[0, 10000], [0,4]])
do_2d_plot(h, (ex, ey), "Signal (ADU)", "Gain Bit Value")
```
%% Cell type:markdown id: tags:
### Mean RAW Preview ###
The per pixel mean of the sequence file of RAW data
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.mean(rim_data,axis=0),
vmin=min(0.75*np.median(rim_data[rim_data > 0]), 2000),
vmax=max(1.5*np.median(rim_data[rim_data > 0]), 16000), cmap="jet")
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:markdown id: tags:
### Mean CORRECTED Preview ###
The per pixel mean of the sequence file of CORR data
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.mean(fim_data,axis=0),
vmin=min(0.75*np.median(fim_data[fim_data > 0]), -0.5),
vmax=max(2.*np.median(fim_data[fim_data > 0]), 100), cmap="jet")
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:markdown id: tags:
### Single Train Preview ###
A single image from the first train
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(fim_data[0,...],
vmin=min(0.75*np.median(fim_data[0,...]), -0.5),
vmax=max(2.*np.median(fim_data[0,...]), 100), cmap="jet")
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:markdown id: tags:
## Signal Distribution ##
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(211)
h = ax.hist(fim_data.flatten(), bins=1000, range=(-100, 1000), log=True)
l = ax.set_xlabel("Signal (keV)")
l = ax.set_ylabel("Counts")
ax = fig.add_subplot(212)
h = ax.hist(fim_data.flatten(), bins=1000, range=(-1000, 10000), log=True)
l = ax.set_xlabel("Signal (keV)")
l = ax.set_ylabel("Counts")
```
%% Cell type:markdown id: tags:
### Maximum GAIN Preview ###
The per pixel maximum of the first train of the GAIN data
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.max(gim_data, axis=0), vmin=0,
vmax=3, cmap="jet")
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:markdown id: tags:
## Bad Pixels ##
The mask contains dedicated entries for all pixels and memory cells as well as all three gains stages. Each mask entry is encoded in 32 bits as:
%% Cell type:code id: tags:
``` python
table = []
for item in BadPixels:
table.append((item.name, "{:016b}".format(item.value)))
md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=["Bad pixel type", "Bit mask"])))
```
%% Cell type:markdown id: tags:
### Single Image Bad Pixels ###
A single image bad pixel map for the first image of the first train
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.log2(msk_data[0,...]), vmin=0, vmax=32, cmap="jet")
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:code id: tags:
``` python
```
......
notebooks/Jungfrau/Jungfrau_dark_analysis_all_gains_burst_mode_NBC.ipynb
+ 60
− 40
View file @ 655e14b1
%% Cell type:markdown id: tags:
# Jungfrau Dark Image Characterization #
Version: 0.1, Author: M. Ramilli, S. Hauf
Analyzes Jungfrau dark image data to deduce offset, noise and resulting bad pixel maps
%% Cell type:code id: tags:
``` python
in_folder = '/gpfs/exfel/exp/FXE/201931/p900088/raw/' # folder under which runs are located, required
cluster_profile = 'noDB' # the ipcluster profile name
in_folder = '/gpfs/exfel/exp/FXE/201931/p900089/raw/' # folder under which runs are located, required
out_folder = '/gpfs/exfel/data/scratch/karnem/test_dark/' # path to place reports at, required
sequences = 1 # number of sequence files in that run
run_high = 86 # run number for G0 dark run, required
run_med = 87 # run number for G1 dark run, required
run_low = 88 # run number for G2 dark run, required
karabo_da = ['JNGFR01'] # list of data aggregators, which corresponds to different JF modules
karabo_id = "FXE_XAD_JF1M" # bla karabo prefix of Jungfrau devices
karabo_id_control = "" # if control is on a different ID, set to empty string if it is the same a karabo-id
receiver_id = 'RECEIVER-{}' # inset for receiver devices
receiver_control_id = "CONTROL" # inset for control devices
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # template to use for file name, double escape sequence number
path_inset = "JNGFR02" # file inset for image data
path_inset_control = "JNGFR02" # file inset for control data
cluster_profile = 'noDB' # the ipcluster profile name
h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data' # path in H5 file under which images are located
h5path_run = '/RUN/{}/DET/{}' # path to run data
h5path_cntrl = '/CONTROL/{}/DET/{}' # path to control data
karabo_da_control = "JNGFR01" # file inset for control data
use_dir_creation_date = True # use dir creation date
cal_db_interface = 'tcp://max-exfl016:8016' # calibrate db interface to connect to
cal_db_timeout = 300000 # timeout on caldb requests
local_output = True # output constants locally
db_output = False # output constants to database
integration_time = 1000 # integration time in us, will be overwritten by value in file
bias_voltage = 90 # sensor bias voltage in V, will be overwritten by value in file
badpixel_threshold_sigma = 20. # bad pixels defined by values outside n times this std from median
offset_abs_threshold_low = [1000, 10000, 10000] # absolute bad pixel threshold in terms of offset, lower values
offset_abs_threshold_high = [8000, 15000, 15000] # absolute bad pixel threshold in terms of offset, upper values
chunkSize = 10 # iteration chunk size, needs to match or be less than number of images in a sequence file
imageRange = [0, 500] # image range in which to evaluate
memoryCells = 16 # number of memory cells
h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data' # path in H5 file under which images are located
run_high = 144 # run number for G0 dark run, required
run_med = 145 # run number for G1 dark run, required
run_low = 146 # run number for G2 dark run, required
karabo_id = "FXE_XAD_JF500K" # karabo prefix of Jungfrau devices
karabo_id_control = "" # if control is on a different ID, set to empty string for using the same as for data
receiver_id = "RECEIVER" # inset for receiver devices
control_id = "CONTROL" # inset for control devices
db_module = "Jungfrau_M260" # ID of module in calibration database
use_dir_creation_date = True # use dir creation date
db_module = ["Jungfrau_M260"] # ID of module in calibration database
manual_slow_data = False # if true, use manually entered bias_voltage and integration_time values
time_limits = 0.025 # to find calibration constants later on, the integration time is allowed to vary by 0.5 us
local_output = True # output constants locally
db_output = False # output constants to database
```
%% Cell type:code id: tags:
``` python
import warnings
warnings.filterwarnings('ignore')
import matplotlib
matplotlib.use('agg')
import numpy as np
import h5py
from h5py import File as h5file
from cal_tools.enums import BadPixels
from cal_tools.tools import get_dir_creation_date, get_random_db_interface
from cal_tools.ana_tools import save_dict_to_hdf5
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from XFELDetAna.util import env
env.iprofile = cluster_profile
from XFELDetAna.detectors.jungfrau import readerPSI as jfreaderPSI
from XFELDetAna.detectors.jungfrau import reader as jfreader
from XFELDetAna.detectors.jungfrau.jf_chunk_reader import JFChunkReader
from XFELDetAna.detectors.jungfrau.util import non_empty_trains, count_n_files, rollout_data, sanitize_data_cellid
import glob
import matplotlib.pyplot as plt
%matplotlib inline
from XFELDetAna.plotting.histogram import histPlot
from XFELDetAna.plotting.heatmap import heatmapPlot
import os
os.makedirs(out_folder, exist_ok=True)
```
%% Cell type:code id: tags:
``` python
path_inset = karabo_da[0] # karabo_da is a concurrency parameter
receiver_id = receiver_id.format(int(path_inset[-2:]))
proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
file_loc = 'proposal:{} runs:{} {} {}'.format(proposal, run_high, run_med, run_low)
# TODO
# this trick is needed until proper mapping is introduced
if len(db_module)>1:
db_module = db_module[int(path_inset[-2:])-1]
else:
db_module = db_module[0]
# Constants relevant for the analysis
run_nums = [run_high, run_med, run_low] # run number for G0/HG0, G1, G2
sensorSize = [1024, 512]
blockSize = [1024, 512]
xRange = [0, 0+sensorSize[0]]
yRange = [0, 0+sensorSize[1]]
gains = [0, 1, 2]
h5path = h5path.format(karabo_id, receiver_id)
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run_high)
print("Using {} as creation time".format(creation_time))
cal_db_interface = get_random_db_interface(cal_db_interface)
print('Calibration database interface: {}'.format(cal_db_interface))
offset_abs_threshold = [offset_abs_threshold_low, offset_abs_threshold_high]
if karabo_id_control == "":
karabo_id_control = karabo_id
import os
os.makedirs(out_folder, exist_ok=True)
```
%% Cell type:code id: tags:
``` python
proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
file_loc = 'proposal:{} runs:{} {} {}'.format(proposal, run_high, run_med, run_low)
print('Path inset ', path_inset)
print('Receiver Id ', receiver_id)
```
%% Cell type:code id: tags:
``` python
def check_memoryCells(file_name, path):
with h5file(file_name, 'r') as f:
t_stamp = np.array(f[path + '/storageCells/timestamp'])
st_cells = np.array(f[path + '/storageCells/value'])
sc_start = np.array(f[path + '/storageCellStart/value'])
valid_train = t_stamp > 0
n_scs = st_cells[valid_train][0] + 1
sc_s = sc_start[valid_train][0]
return n_scs, sc_s
```
%% Cell type:code id: tags:
``` python
chunkSize = 100
filep_size = 1000
noiseCal = None
noise_map = None
offset_map = None
memoryCells = None
for i, r_n in enumerate(run_nums):
gain = i
print(f"Gain stage {gain}, run {r_n}")
valid_data = []
valid_cellids = []
if r_n is not None:
n_tr = 0
n_empty_trains = 0
n_empty_sc = 0
ped_dir = "{}/r{:04d}/".format(in_folder, r_n)
fp_name = path_template.format(r_n, path_inset_control)
fp_name = path_template.format(r_n, karabo_da_control)
fp_path = '{}/{}'.format(ped_dir, fp_name)
n_files = len(glob.glob("{}/*{}*.h5".format(ped_dir, path_inset)))
myRange = range(0, n_files)
control_path = '/CONTROL/{}/DET/{}'.format(karabo_id_control, control_id)
control_path = h5path_cntrl.format(karabo_id_control, receiver_control_id)
this_run_mcells, sc_start = check_memoryCells(fp_path.format(0).format(myRange[0]), control_path)
if noise_map is None:
if not manual_slow_data:
with h5py.File(fp_path.format(0), 'r') as f:
integration_time = float(f['/RUN/{}/DET/{}/exposureTime/value'.format(karabo_id_control, control_id)][()]*1e6)
bias_voltage = int(np.squeeze(f['/RUN/{}/DET/{}/vHighVoltage/value'.format(karabo_id_control, control_id)])[0])
run_path = h5path_run.format(karabo_id_control, receiver_control_id)
integration_time = float(f[f'{run_path}/exposureTime/value'][()]*1e6)
bias_voltage = int(np.squeeze(f[f'{run_path}/vHighVoltage/value'])[0])
print("Integration time is {} us".format(integration_time))
print("Bias voltage is {} V".format(bias_voltage))
if this_run_mcells == 1:
memoryCells = 1
print('Dark runs in single cell mode\n storage cell start: {:02d}'.format(sc_start))
else:
memoryCells = 16
print('Dark runs in burst mode\n storage cell start: {:02d}'.format(sc_start))
noise_map = np.zeros(sensorSize+[memoryCells, 3])
offset_map = np.zeros(sensorSize+[memoryCells, 3])
fp_name = path_template.format(r_n, path_inset)
fp_path = '{}/{}'.format(ped_dir, fp_name)
myRange_P = range(0, sequences)
path = h5path
print("Reading data from {}".format(fp_path))
print("Run is: {}".format(r_n))
print("HDF5 path: {}".format(h5path))
imageRange = [0, filep_size*len(myRange)]
reader = JFChunkReader(filename = fp_path, readFun = jfreader.readData, size = filep_size, chunkSize = chunkSize,
path = h5path, image_range=imageRange, pixels_x = sensorSize[0], pixels_y = sensorSize[1],
x_range = xRange, y_range = yRange, imagesPerChunk=chunkSize, filesRange = myRange,
memoryCells=this_run_mcells, blockSize=blockSize)
for data in reader.readChunks():
images = np.array(data[0], dtype=np.float)
gainmaps = np.array(data[1], dtype=np.uint16)
trainId = np.array(data[2])
fr_num = np.array(data[3])
acelltable = np.array(data[4])
n_tr += acelltable.shape[-1]
this_tr = acelltable.shape[-1]
idxs = non_empty_trains(trainId)
images = images[..., idxs]
gainmaps = gainmaps[..., idxs]
fr_num = fr_num[..., idxs]
acelltable = acelltable[..., idxs]
if memoryCells == 1:
acelltable -= sc_start
n_empty_trains += this_tr - acelltable.shape[-1]
n_empty_sc += len(acelltable[acelltable > 15])
if i > 0 and memoryCells == 16: ## throwing away the second memory cell entry in lower gains
acelltable[1] = 255
images, gainmaps, acelltable = rollout_data([images, gainmaps, acelltable]) # makes 4-dim vecs into 3-dim
# makes 2-dim into 1-dim
# leaves 1-dim and 3-dim vecs
images, gainmaps, acelltable = sanitize_data_cellid([images, gainmaps], acelltable) # removes entries with cellID 255
valid_data.append(images)
valid_cellids.append(acelltable)
valid_data = np.concatenate(valid_data, axis=2)
valid_cellids = np.concatenate(valid_cellids, axis=0)
for cell in range(memoryCells):
thiscell = valid_data[...,valid_cellids == cell]
noise_map[...,cell,gain] = np.std(thiscell, axis=2)
offset_map[...,cell,gain] = np.mean(thiscell, axis=2)
print('G{:01d} dark calibration'.format(i))
print('Missed {:d} out of {:d} trains'.format(n_empty_trains, n_tr))
print('Lost {:d} images out of {:d}'.format(n_empty_sc, this_run_mcells * (n_tr - n_empty_trains)))
else:
print('missing G{:01d}'.format(i))
```
%% Cell type:markdown id: tags:
## Offset and Noise Maps ##
Below offset and noise maps for the high ($g_0$) gain stage are shown, alongside the distribution of these values. One expects block-like structures mapping to the ASICs of the detector
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
from XFELDetAna.core.util import remove_nans
%matplotlib inline
#%matplotlib notebook
from XFELDetAna.plotting.histogram import histPlot
from XFELDetAna.plotting.heatmap import heatmapPlot
g_name = ['G0', 'G1', 'G2']
g_range = [(0, 8000), (8000, 16000), (8000, 16000)]
n_range = [(0., 50.), (0., 50.), (0., 50.)]
unit = '[ADCu]'
```
%% Cell type:code id: tags:
``` python
for g_idx in gains:
for cell in range(0, memoryCells):
f_o0 = heatmapPlot(np.swapaxes(offset_map[..., cell, g_idx], 0, 1),
y_label="Row",
x_label="Column",
lut_label="Pedestal {} [ADCu]".format(g_name[g_idx]),
aspect=1.,
vmin=g_range[g_idx][0],
vmax=g_range[g_idx][1])
fo0, ax_o0 = plt.subplots()
res_o0 = histPlot(ax_o0, offset_map[..., cell, g_idx],
bins=800,
range=g_range[g_idx],
facecolor='b',
histotype='stepfilled')
ax_o0.tick_params(axis='both',which='major',labelsize=15)
ax_o0.set_title('Module pedestal distribution Cell {:02d}'.format(cell), fontsize=15)
ax_o0.set_xlabel('Pedestal {} [ADCu]'.format(g_name[g_idx]),fontsize=15)
ax_o0.set_yscale('log')
f_n0 = heatmapPlot(np.swapaxes(noise_map[..., cell, g_idx], 0, 1),
y_label="Row",
x_label="Column",
lut_label="RMS noise {} ".format(g_name[g_idx]) + unit,
aspect=1.,
vmin=n_range[g_idx][0],
vmax=n_range[g_idx][1])
fn0, ax_n0 = plt.subplots()
res_n0 = histPlot(ax_n0, noise_map[..., cell, g_idx],
bins=100,
range=n_range[g_idx],
facecolor='b',
histotype='stepfilled')
ax_n0.tick_params(axis='both',which='major',labelsize=15)
ax_n0.set_title('Module noise distribution Cell {:02d}'.format(cell), fontsize=15)
ax_n0.set_xlabel('RMS noise {} '.format(g_name[g_idx]) + unit,fontsize=15)
#ax_n0.set_yscale('log')
plt.show()
```
%% Cell type:markdown id: tags:
## Bad Pixel Map ###
The bad pixel map is deduced by comparing offset and noise of each pixel ($v_i$) and each gain ($g$) against the median value for that gain stage:
$$
v_i > \mathrm{median}(v_{k,g}) + n \sigma_{v_{k,g}}
$$
or
$$
v_i < \mathrm{median}(v_{k,g}) - n \sigma_{v_{k,g}}
$$
Values are encoded in a 32 bit mask, where for the dark image deduced bad pixels the following non-zero entries are relevant:
%% Cell type:code id: tags:
``` python
def print_bp_entry(bp):
print("{:<30s} {:032b}".format(bp.name, bp.value))
print_bp_entry(BadPixels.OFFSET_OUT_OF_THRESHOLD)
print_bp_entry(BadPixels.NOISE_OUT_OF_THRESHOLD)
print_bp_entry(BadPixels.OFFSET_NOISE_EVAL_ERROR)
```
%% Cell type:code id: tags:
``` python
bad_pixels_map = np.zeros(noise_map.shape, np.uint32)
def eval_bpidx(d):
mdn = np.nanmedian(d, axis=(0, 1, 2))[None, None, None, :]
std = np.nanstd(d, axis=(0, 1, 2))[None, None, None, :]
idx = (d > badpixel_threshold_sigma*std+mdn) | (d < (-badpixel_threshold_sigma)*std+mdn)
return idx
offset_abs_threshold = np.array(offset_abs_threshold)
bad_pixels_map[eval_bpidx(offset_map)] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value
bad_pixels_map[~np.isfinite(offset_map)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
bad_pixels_map[eval_bpidx(noise_map)] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
bad_pixels_map[~np.isfinite(noise_map)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
bad_pixels_map[(offset_map < offset_abs_threshold[0][None, None, None, :]) | (offset_map > offset_abs_threshold[1][None, None, None, :])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value
for g_idx in gains:
for cell in range(0, memoryCells):
bad_pixels = bad_pixels_map[:, :, cell, g_idx]
fn_0 = heatmapPlot(np.swapaxes(bad_pixels, 0, 1),
y_label="Row",
x_label="Column",
lut_label="Badpixels {} [ADCu]".format(g_name[g_idx]),
aspect=1.,
vmin=0)
```
%% Cell type:code id: tags:
``` python
constants = {'Offset': np.moveaxis(offset_map, 0, 1),
'Noise': np.moveaxis(noise_map, 0, 1),
'BadPixelsDark': np.moveaxis(bad_pixels_map, 0, 1)}
for key, const_data in constants.items():
metadata = ConstantMetaData()
const = getattr(Constants.jungfrau, key)()
const.data = const_data
metadata.calibration_constant = const
# set the operating condition
condition = Conditions.Dark.jungfrau(memory_cells=memoryCells, bias_voltage=bias_voltage,
integration_time=integration_time)
for parm in condition.parameters:
if parm.name == "Integration Time":
parm.lower_deviation = time_limits
parm.upper_deviation = time_limits
device = getattr(Detectors, db_module)
metadata.detector_condition = condition
# specify the 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.calibration_constant_version.raw_data_location = file_loc
if db_output:
metadata.send(cal_db_interface)
metadata.send(cal_db_interface, timeout=cal_db_timeout)
print('Constants {} is sent to the data base'.format(key))
# save everything to file.
# one file - one entrt in the DB
if local_output:
dpar = {}
for parm in metadata.detector_condition.parameters:
dpar[parm.name] = {'lower_deviation': parm.lower_deviation,
'upper_deviation': parm.upper_deviation,
'value': parm.value}
data_to_store = {}
data_to_store['condition'] = dpar
data_to_store['db_module'] = db_module
data_to_store['constant'] = key
data_to_store['data'] = const_data
data_to_store['creation_time'] = creation_time
data_to_store['dclass'] = 'jungfrau'
data_to_store['file_loc'] = file_loc
ofile = "{}/const_{}_{}.h5".format(out_folder, key, db_module)
save_dict_to_hdf5(data_to_store, ofile)
```
......
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