Skip to content
Snippets Groups Projects
Commit b8eafa40 authored by Cyril Danilevski's avatar Cyril Danilevski :scooter:
Browse files

Merge branch 'feat/jungfrau_gain_setting' into 'master' 

[Jungfrau] Add Gain setting to Jungfrau notebooks

See merge request detectors/pycalibration!555
parents 777c5c55 5c857454
No related branches found
No related tags found
1 merge request!555[Jungfrau] Add Gain setting to Jungfrau notebooks
Hide whitespace changes
Inline Side-by-side
notebooks/Jungfrau/Jungfrau_Gain_Correct_and_Verify_NBC.ipynb
+ 7
− 3
View file @ b8eafa40
%% 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/CALLAB/202031/p900113/raw" # the folder to read data from, required
out_folder = "/gpfs/exfel/data/scratch/hammerd/issue-242" # the folder to output to, required
sequences = [-1] # sequences to correct, set to [-1] for all, range allowed
run = 9979 # run to process, required
karabo_id = "SPB_IRDA_JF4M" # karabo prefix of Jungfrau devices
karabo_da = ['JNGFR01'] # data aggregators
receiver_id = "JNGFR{:02d}" # 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
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 = "JNGFRCTRL00" # 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-exfl016:8017#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
bias_voltage = 180 # will be overwritten by value in file
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.
integration_time = 4.96 # integration time in us, will be overwritten by value in file
gain_setting = 0 # 0 for dynamic gain, 1 for dynamic HG0, will be overwritten by value in file
mem_cells = 0 # leave memory cells equal 0, as it is saved in control information starting 2019.
db_module = ["Jungfrau_M275"] # 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, karabo_da):
from xfel_calibrate.calibrate import balance_sequences as bs
return bs(in_folder, run, sequences, sequences_per_node, karabo_da)
```
%% Cell type:code id: tags:
``` python
import copy
import multiprocessing
import time
import warnings
from functools import partial
from pathlib import Path
import h5py
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import tabulate
from cal_tools.enums import BadPixels
from cal_tools.tools import (
get_constant_from_db_and_time,
get_dir_creation_date,
map_modules_from_folder,
)
from iCalibrationDB import Conditions, Constants
from IPython.display import Latex, display
from matplotlib.colors import LogNorm
warnings.filterwarnings('ignore')
matplotlib.use('agg')
%matplotlib inline
```
%% Cell type:code id: tags:
``` python
in_folder = Path(in_folder)
out_folder = Path(out_folder)
ped_dir = in_folder / f'r{run:04d}'
h5path = h5path.format(karabo_id, receiver_id)
if out_folder.exists() and not overwrite:
raise AttributeError("Output path exists! Exiting")
else:
out_folder.mkdir(parents=True, exist_ok=True)
fp_name_contr = path_template.format(run, karabo_da_control, 0)
fp_path_contr = ped_dir / fp_name_contr
if sequences[0] == -1:
sequences = None
print(f"Run is: {run}")
print(f"HDF5 path: {h5path}")
print(f"Process modules: {karabo_da}")
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run)
print(f"Using {creation_time} as creation time")
if karabo_id_control == "":
karabo_id_control = karabo_id
```
%% Cell type:code id: tags:
``` python
def check_memory_cells(file_name, path):
with h5py.File(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
# set everything up filewise
mapped_files, mod_ids, total_sequences, sequences_qm, _ = map_modules_from_folder(
in_folder, run, path_template, karabo_da, sequences
)
print(f"Processing a total of {total_sequences} sequence files")
table = []
fi = 0
if total_sequences > 0: # create table
for i, key in enumerate(mapped_files):
for k, f in enumerate(list(mapped_files[key].queue)):
if k == 0:
table.append((fi, karabo_da[i], k, f))
else:
table.append((fi, "", k, f))
fi += 1
md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
headers=["#", "module", "# module", "file"])))
# restore the queue
mapped_files, mod_ids, total_sequences, sequences_qm, _ = map_modules_from_folder(
in_folder, run, path_template, karabo_da, sequences
)
```
%% Cell type:code id: tags:
``` python
if not manual_slow_data:
with h5py.File(fp_path_contr, 'r') as f:
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])
gain_s = f[f'/RUN/{karabo_id_control}/DET/CONTROL/settings/value'][0].decode()
gain_setting = 1 if gain_s == 'dynamichg0' else 0 # gain_s == 'dynamicgain'
control_path = h5path_cntrl.format(karabo_id_control, receiver_control_id)
try:
this_run_mcells, sc_start = check_memory_cells(fp_path_contr, control_path)
if this_run_mcells == 1:
memory_cells = 1
print(f'Dark runs in single cell mode\n storage cell start: {sc_start:02d}')
else:
memory_cells = 16
print(f'Dark runs in burst mode\n storage cell start: {sc_start:02d}')
except Exception as e:
if "Unable to open object" in str(e):
if mem_cells==0:
memory_cells = 1
else:
memory_cells = mem_cells
print(f'Set memory cells to {memory_cells} as it is not saved in control information.')
else:
print(f"Error trying to access memory cell from contol information: {e}")
print(f"Integration time is {integration_time} us")
print(f"Gain setting is {gain_setting}")
print(f"Bias voltage is {bias_voltage} V")
print(f"Number of memory cells is {memory_cells}")
```
%% Cell type:code id: tags:
``` python
condition = Conditions.Dark.jungfrau(
memory_cells=memory_cells,
bias_voltage=bias_voltage,
integration_time=integration_time,
gain_setting=gain_setting,
)
def get_constants_for_module(karabo_da: str):
""" Get calibration constants for given module of Jungfrau
:return:
offset_map (offset map),
mask (mask of bad pixels),
gain_map (map of relative gain factors),
db_module (name of DB module),
when (dictionaty: constant - creation time)
"""
when = {}
retrieval_function = partial(
get_constant_from_db_and_time,
karabo_id=karabo_id,
karabo_da=karabo_da,
condition=condition,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout,
)
offset_map, when["Offset"] = retrieval_function(
constant=Constants.jungfrau.Offset(), empty_constant=np.zeros((1024, 512, 1, 3))
)
mask, when["BadPixelsDark"] = retrieval_function(
constant=Constants.jungfrau.BadPixelsDark(),
empty_constant=np.zeros((1024, 512, 1, 3)),
)
mask_ff, when["BadPixelsFF"] = retrieval_function(
constant=Constants.jungfrau.BadPixelsFF(),
empty_constant=None
)
gain_map, when["Gain"] = retrieval_function(
constant=Constants.jungfrau.RelativeGain(),
empty_constant=None
)
# combine masks
if mask_ff is not None:
mask |= np.moveaxis(mask_ff, 0, 1)
# move from x,y,cell,gain to cell,x,y,gain
offset_map = np.squeeze(offset_map)
mask = np.squeeze(mask)
if memory_cells > 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 not None:
if memory_cells > 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)
return offset_map, mask, gain_map, karabo_da, when
with multiprocessing.Pool() as pool:
r = pool.map(get_constants_for_module, karabo_da)
constants = {}
for offset_map, mask, gain_map, k_da, when in r:
print(f'Constants for module {k_da}:')
for const in when:
print(f' {const} injected at {when[const]}')
if gain_map is None:
print(" No gain map found")
no_relative_gain = True
constants[k_da] = (offset_map, mask, gain_map)
```
%% 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`
"""
"""Copy and sanitize data from `infile` that is not calibrated."""
h5base = h5base.lstrip("/")
dont_copy = ["adc", ]
dont_copy = [f'{h5base}/{dnc}' for dnc 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
# Correct a chunk of images for offset and gain
def correct_chunk(offset_map, mask, gain_map, memory_cells, no_relative_gain, inp):
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 memory_cells==1:
rim_data = np.squeeze(copy.copy(d))
else:
rim_data = np.squeeze(copy.copy(d[:,0,...]))
# Select memory cells
if memory_cells>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 memory_cells>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 memory_cells==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,...]))
except Exception as e:
err = e
return ind, d, msk, rim_data, fim_data, gim_data, msk_data, err
```
%% Cell type:code id: tags:
``` python
fim_data = {}
gim_data = {}
rim_data = {}
msk_data = {}
# For each module, chunks will be processed by pool
pool = multiprocessing.Pool()
# Loop over modules
for local_karabo_da, mapped_files_module in zip(karabo_da, mapped_files.values()):
h5path_f = h5path.format(int(local_karabo_da[-2:]))
# Loop over sequences for given module
for sequence_file_number, sequence_file in enumerate(mapped_files_module.queue):
sequence_file = Path(sequence_file)
offset_map, mask, gain_map = constants[local_karabo_da]
with h5py.File(sequence_file, 'r') as infile:
# The processed files are saved here in a folder with the run name.
out_filename = out_folder / sequence_file.name.replace("RAW", "CORR")
print(f'Process file: {sequence_file}, with path {h5path_f}')
try:
with h5py.File(out_filename, "w") as outfile:
copy_and_sanitize_non_cal_data(infile, outfile, h5path_f)
oshape = infile[h5path_f+"/adc"].shape
print(f'Data shape: {oshape}')
if not oshape[0]:
raise ValueError(f"No image data: shape {oshape}")
# Chunk always contains >= 1 complete image
chunk_shape = (chunk_size_idim, 1) + oshape[-2:]
ddset = outfile.create_dataset(h5path_f+"/adc",
oshape,
chunks=chunk_shape,
dtype=np.float32)
mskset = outfile.create_dataset(h5path_f+"/mask",
oshape,
chunks=chunk_shape,
dtype=np.uint32,
compression="gzip", compression_opts=1, shuffle=True)
# Parallelize 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(f'Chunk size: {chunk_size}')
ts = time.time()
while ind<max_ind:
d = infile[h5path_f+"/adc"][ind:ind+chunk_size,...].astype(np.float32)
g = infile[h5path_f+"/gain"][ind:ind+chunk_size,...]
if h5path_f+"/memoryCell" in infile:
m = infile[h5path_f+"/memoryCell"][ind:ind+chunk_size,...]
else:
m = None
print(f'To process: {d.shape}')
inp.append((d, g, m, ind, sequence_file_number==0))
ind += chunk_size
print('Preparation time: ', time.time() - ts)
ts = time.time()
print(f'Run {len(inp)} processes')
p = partial(correct_chunk, offset_map, mask, gain_map, memory_cells, no_relative_gain)
r = pool.map(p, inp)
if sequence_file_number == 0:
(_,_,_,
rim_data[local_karabo_da], fim_data[local_karabo_da],
gim_data[local_karabo_da], msk_data[local_karabo_da], _) = r[0]
print('Correction time: ', time.time() - ts)
ts = time.time()
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(f'Error: {err}')
print('Saving time: ', time.time() - ts)
except Exception as e:
print(f"Error: {e}")
pool.close()
```
%% Cell type:code id: tags:
``` python
def do_2d_plot(data, edges, y_axis, x_axis, title):
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)
ax.set_title(title)
cb = fig.colorbar(im)
cb.set_label("Counts")
```
%% Cell type:code id: tags:
``` python
for mod in rim_data:
h, ex, ey = np.histogram2d(rim_data[mod].flatten(),
gim_data[mod].flatten(),
bins=[100, 4],
range=[[0, 10000], [0, 4]])
do_2d_plot(h, (ex, ey), "Signal (ADU)", "Gain Bit Value", f'Module {mod}')
```
%% 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
for mod in rim_data:
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.mean(rim_data[mod],axis=0),
vmin=min(0.75*np.median(rim_data[mod][rim_data[mod] > 0]), 2000),
vmax=max(1.5*np.median(rim_data[mod][rim_data[mod] > 0]), 16000), cmap="jet")
ax.set_title(f'Module {mod}')
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
for mod in rim_data:
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.mean(fim_data[mod],axis=0),
vmin=min(0.75*np.median(fim_data[mod][fim_data[mod] > 0]), -0.5),
vmax=max(2.*np.median(fim_data[mod][fim_data[mod] > 0]), 100), cmap="jet")
ax.set_title(f'Module {mod}', size=18)
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
for mod in rim_data:
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(fim_data[mod][0,...],
vmin=min(0.75*np.median(fim_data[mod][0,...]), -0.5),
vmax=max(2.*np.median(fim_data[mod][0,...]), 100), cmap="jet")
ax.set_title(f'Module {mod}', size=18)
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:markdown id: tags:
## Signal Distribution ##
%% Cell type:code id: tags:
``` python
for mod in rim_data:
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(211)
h = ax.hist(fim_data[mod].flatten(), bins=1000, range=(-100, 1000), log=True)
l = ax.set_xlabel("Signal (keV)")
l = ax.set_ylabel("Counts")
_ = ax.set_title(f'Module {mod}')
ax = fig.add_subplot(212)
h = ax.hist(fim_data[mod].flatten(), bins=1000, range=(-1000, 10000), log=True)
l = ax.set_xlabel("Signal (keV)")
l = ax.set_ylabel("Counts")
_ = ax.set_title(f'Module {mod}')
```
%% 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
for mod in rim_data:
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.max(gim_data[mod], axis=0), vmin=0,
vmax=3, cmap="jet")
ax.set_title(f'Module {mod}', size=18)
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, f"{item.value:016b}"))
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
for mod in rim_data:
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.log2(msk_data[mod][0,...]), vmin=0, vmax=32, cmap="jet")
ax.set_title(f'Module {mod}', size=18)
cb = fig.colorbar(im, ax=ax)
```
......
notebooks/Jungfrau/Jungfrau_dark_analysis_all_gains_burst_mode_NBC.ipynb
+ 12
− 4
View file @ b8eafa40
%% 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
cluster_profile = 'noDB' # the ipcluster profile name
in_folder = '/gpfs/exfel/exp/SPB/202130/p900204/raw/' # folder under which runs are located, required
out_folder = '/gpfs/exfel/data/scratch/ahmedk/jftest_dark/' # path to place reports at, required
run_high = 141 # run number for G0 dark run, required
run_med = 142 # run number for G1 dark run, required
run_low = 143 # run number for G2 dark run, required
karabo_da = ['JNGFR01', 'JNGFR02','JNGFR03','JNGFR04', 'JNGFR05', 'JNGFR06','JNGFR07','JNGFR08'] # list of data aggregators, which corresponds to different JF modules
karabo_id = "SPB_IRDA_JF4M" # karabo_id (detector identifier) prefix of Jungfrau detector to process.
karabo_id_control = "" # if control is on a different ID, set to empty string if it is the same a karabo-id
receiver_id = 'JNGFR{:02}' # 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_da_control = "JNGFRCTRL00" # 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
gain_setting = 0 # 0 for dynamic, forceswitchg1, forceswitchg2, 1 for dynamichg0, fixedgain1, fixgain2. 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 = 5. # 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
db_module = ['Jungfrau_M275', "Jungfrau_M035", 'Jungfrau_M273','Jungfrau_M203','Jungfrau_M221','Jungfrau_M267'] # 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
operation_mode = '' # Detector operation mode, optional
```
%% Cell type:code id: tags:
``` python
import glob
import os
import warnings
warnings.filterwarnings('ignore')
import h5py
import matplotlib
from h5py import File as h5file
matplotlib.use('agg')
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
from cal_tools.ana_tools import save_dict_to_hdf5
from cal_tools.enums import BadPixels
from cal_tools.tools import (
get_dir_creation_date,
get_pdu_from_db,
get_random_db_interface,
get_report,
save_const_to_h5,
send_to_db,
)
from iCalibrationDB import Conditions, Constants, Detectors, Versions
from XFELDetAna.util import env
env.iprofile = cluster_profile
from XFELDetAna.detectors.jungfrau import reader as jfreader
from XFELDetAna.detectors.jungfrau import readerPSI as jfreaderPSI
from XFELDetAna.detectors.jungfrau.jf_chunk_reader import JFChunkReader
from XFELDetAna.detectors.jungfrau.util import (
count_n_files,
rollout_data,
sanitize_data_cellid,
)
from XFELDetAna.plotting.heatmap import heatmapPlot
from XFELDetAna.plotting.histogram import histPlot
```
%% 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)
report = get_report(out_folder)
os.makedirs(out_folder, exist_ok=True)
# TODO
# this trick is needed until proper mapping is introduced
if len(db_module)>1:
# TODO: SPB JF Hack till using all modules.
if karabo_id == "SPB_IRDA_JNGFR" and int(path_inset[-2:]) > 5:
db_module = db_module[int(path_inset[-2:])-3]
else:
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
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, karabo_da_control)
fp_path = '{}/{}'.format(ped_dir, fp_name)
files_pattern = "{}/*{}*.h5".format(ped_dir, path_inset)
n_files = len(glob.glob(files_pattern))
if n_files == 0:
raise Exception(f"No files found matching {files_pattern!r}")
myRange = range(0, n_files)
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:
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])
if r_n == run_high:
gain_s = f[f'/RUN/{karabo_id_control}/DET/CONTROL/settings/value'][0].decode()
gain_setting = 1 if gain_s == "dynamichg0" else 0
print(f"Constants Gain setting is {gain_setting} ({gain_s})")
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)
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 = np.nonzero(trainId)[0]
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 all the SC entries except the first for lower gains
acelltable[1:] = 255
# makes 4-dim vecs into 3-dim
# makes 2-dim into 1-dim
# leaves 1-dim and 3-dim vecs
images, gainmaps, acelltable = rollout_data([images, gainmaps, acelltable])
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.heatmap import heatmapPlot
from XFELDetAna.plotting.histogram import histPlot
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=unit,
aspect=1.,
vmin=g_range[g_idx][0],
vmax=g_range[g_idx][1],
title=f'Pedestal {g_name[g_idx]} - Cell {cell:02d}')
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(f'Module pedestal distribution - Cell {cell:02d}', fontsize=15)
ax_o0.set_xlabel(f'Pedestal {g_name[g_idx]} {unit}',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= unit,
aspect=1.,
vmin=n_range[g_idx][0],
vmax=n_range[g_idx][1],
title=f"RMS noise {g_name[g_idx]} - Cell {cell:02d}")
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(f'Module noise distribution - Cell {cell:02d}', fontsize=15)
ax_n0.set_xlabel(f'RMS noise {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))[None, None, :, :]
std = np.nanstd(d, axis=(0, 1))[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(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=f"Badpixels {g_name[g_idx]} [ADCu]",
aspect=1.,
vmin=0, title=f'G{g_idx} Bad pixel map - Cell {cell:02d}')
```
%% Cell type:code id: tags:
``` python
# TODO: this cell in the notebook is not designed to run for more than one module
# Constants need to be able to have constant for each module as for big detectors
constants = {'Offset': np.moveaxis(offset_map, 0, 1),
'Noise': np.moveaxis(noise_map, 0, 1),
'BadPixelsDark': np.moveaxis(bad_pixels_map, 0, 1)}
md = None
for key, const_data in constants.items():
const = getattr(Constants.jungfrau, key)()
const.data = const_data
# set the operating condition
condition = Conditions.Dark.jungfrau(memory_cells=memoryCells, bias_voltage=bias_voltage,
integration_time=integration_time)
integration_time=integration_time,
gain_setting=gain_setting)
for parm in condition.parameters:
if parm.name == "Integration Time":
parm.lower_deviation = time_limits
parm.upper_deviation = time_limits
# This should be used in case of running notebook
# by a different method other than myMDC which already
# sends CalCat info.
# TODO: Set db_module to "" by default in the first cell
if not db_module:
db_module = get_pdu_from_db(karabo_id, karabo_da, const,
condition, cal_db_interface,
snapshot_at=creation_time)[0]
if db_output:
md = send_to_db(db_module, karabo_id, const, condition,
file_loc=file_loc, report_path=report,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout)
if local_output:
md = save_const_to_h5(db_module, karabo_id, const, condition,
const.data, file_loc, report,
creation_time, out_folder)
print(f"Calibration constant {key} is stored locally at {out_folder}.\n")
print("Constants parameter conditions are:\n")
print(f"• Bias voltage: {bias_voltage}\n• Memory cells: {memoryCells}\n"
f"• Integration time: {integration_time}\n"
f"• Gain setting: {gain_setting}\n"
f"• Creation time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\n")
```
%% Cell type:code id: tags:
``` python
```
......
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