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

Merge branch 'deployment_issues' into 'master' 

Fix issues found in testing

See merge request detectors/pycalibration!584
parents 11a201bc 5c208f82
No related branches found
No related tags found
1 merge request!584Fix issues found in testing
Hide whitespace changes
Inline Side-by-side
notebooks/Jungfrau/Jungfrau_dark_analysis_all_gains_burst_mode_NBC.ipynb
+ 3
− 2
View file @ f77e1e5d
%% 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 = '' # 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 = "" # ID of module in calibration database, this parameter is ignored in the notebook. TODO: remove from calibration_configurations.
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
import matplotlib.pyplot as plt
import numpy as np
from h5py import File as h5file
matplotlib.use('agg')
%matplotlib inline
from XFELDetAna.detectors.jungfrau.util import (
rollout_data,
sanitize_data_cellid,
)
from XFELDetAna.plotting.heatmap import heatmapPlot
from XFELDetAna.plotting.histogram import histPlot
from XFELDetAna.detectors.jungfrau import reader as jfreader
from XFELDetAna.detectors.jungfrau.jf_chunk_reader import JFChunkReader
from XFELDetAna.util import env
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
env.iprofile = cluster_profile
```
%% 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)
# 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
memoryCells = None
noise_map = dict()
offset_map = dict()
# TODO: parallelize with multiprocessing.
for mod in karabo_da:
for gain, r_n in enumerate(run_nums):
print(f"Gain stage {gain}, run {r_n}")
valid_data = []
valid_cellids = []
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 mod not in noise_map:
if not manual_slow_data:
run_path = h5path_run.format(karabo_id_control, receiver_control_id)
with h5py.File(fp_path.format(0), 'r') as f:
filename = fp_path.format(0)
with h5py.File(filename, 'r') as f:
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:
try:
gain_s = f[f'/RUN/{karabo_id_control}/DET/CONTROL/settings/value'][0].decode()
except KeyError:
print(
"ERROR: gain_setting is not available for h5 ctrl path "
f"/RUN/{karabo_id_control}/DET/CONTROL/settings/value,\nfor file: {fp_path}. \n"
f"/RUN/{karabo_id_control}/DET/CONTROL/settings/value,\nfor file: {filename}. \n"
"WARNING: Setting gain_setting to 0, assuming that this is an old run.\n")
gain_s = "KeyError"
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[mod] = np.zeros(sensorSize+[memoryCells, 3])
offset_map[mod] = np.zeros(sensorSize+[memoryCells, 3])
fp_name = path_template.format(r_n, path_inset)
fp_path = '{}/{}'.format(ped_dir, fp_name)
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])
# throwing away all the SC entries except
# the first for lower gains.
if gain > 0 and memoryCells == 16:
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])
# removes entries with cellID 255
images, gainmaps, acelltable = sanitize_data_cellid([images, gainmaps], acelltable)
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[mod][..., cell, gain] = np.std(thiscell, axis=2)
offset_map[mod][..., cell, gain] = np.mean(thiscell, axis=2)
print(f'G{gain:01d} dark calibration')
print(f'Missed {n_empty_trains:d} out of {n_tr:d} trains')
print(
f'Lost {n_empty_sc:d} images out of {this_run_mcells*(n_tr-n_empty_trains):d}') # noqa
```
%% 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
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 mod in karabo_da:
for g_idx in gains:
for cell in range(0, memoryCells):
f_o0 = heatmapPlot(
np.swapaxes(offset_map[mod][..., 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} - Module {mod}')
fo0, ax_o0 = plt.subplots()
res_o0 = histPlot(
ax_o0, offset_map[mod][..., 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} - Module {mod}',
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[mod][..., 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} - Module {mod}",
)
fn0, ax_n0 = plt.subplots()
res_n0 = histPlot(
ax_n0,
noise_map[mod][..., 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} - Module {mod}',
fontsize=15)
ax_n0.set_xlabel(
f'RMS noise {g_name[g_idx]} ' + unit, fontsize=15)
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
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
```
%% Cell type:code id: tags:
``` python
bad_pixels_map = dict()
for mod in karabo_da:
bad_pixels_map[mod] = np.zeros(noise_map[mod].shape, np.uint32)
offset_abs_threshold = np.array(offset_abs_threshold)
bad_pixels_map[mod][eval_bpidx(offset_map[mod])] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value
bad_pixels_map[mod][~np.isfinite(offset_map[mod])] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
bad_pixels_map[mod][eval_bpidx(noise_map[mod])] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
bad_pixels_map[mod][~np.isfinite(noise_map[mod])] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
bad_pixels_map[mod][(offset_map[mod] < offset_abs_threshold[0][None, None, None, :]) | (offset_map[mod] > offset_abs_threshold[1][None, None, None, :])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value # noqa
for g_idx in gains:
for cell in range(memoryCells):
bad_pixels = bad_pixels_map[mod][:, :, 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} - Module {mod}')
```
%% Cell type:code id: tags:
``` python
# set the operating condition
condition = Conditions.Dark.jungfrau(
memory_cells=memoryCells,
bias_voltage=bias_voltage,
integration_time=integration_time,
gain_setting=gain_setting)
db_modules = get_pdu_from_db(
karabo_id=karabo_id,
karabo_da=karabo_da,
constant=Constants.jungfrau.Offset(),
condition=condition,
cal_db_interface=cal_db_interface,
snapshot_at=creation_time)
```
%% Cell type:code id: tags:
``` python
for mod, db_mod in zip(karabo_da, db_modules):
constants = {
'Offset': np.moveaxis(offset_map[mod], 0, 1),
'Noise': np.moveaxis(noise_map[mod], 0, 1),
'BadPixelsDark': np.moveaxis(bad_pixels_map[mod], 0, 1),
}
md = None
for key, const_data in constants.items():
const = getattr(Constants.jungfrau, key)()
const.data = const_data
for parm in condition.parameters:
if parm.name == "Integration Time":
parm.lower_deviation = time_limits
parm.upper_deviation = time_limits
if db_output:
md = send_to_db(
db_module=db_mod,
karabo_id=karabo_id,
constant=const,
condition=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=db_mod,
karabo_id=karabo_id,
constant=const,
condition=condition,
data=const.data,
file_loc=file_loc,
report=report,
creation_time=creation_time,
out_folder=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"
f"• 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") # noqa
```
......
src/xfel_calibrate/calibrate.py
+ 1
− 1
View file @ f77e1e5d
......@@ -717,7 +717,7 @@ class SlurmOptions:
launcher_slurm += ["--job-name", self.job_name]
if self.nice:
launcher_slurm += ["--nice", str(self.nice)]
launcher_slurm.append(f"--nice={self.nice}")
if self.mem:
launcher_slurm.append(f"--mem={self.mem}G")
......
webservice/webservice.py
+ 16
− 7
View file @ f77e1e5d
......@@ -513,13 +513,22 @@ async def run_action(job_db, cmd, mode, proposal, run, rid) -> str:
rstr = stdout.decode()
for r in rstr.split("\n"):
if "Submitted job:" in r:
_, jobid = r.split(":")
c.execute(
"INSERT INTO jobs VALUES (?, ?, ?, ?, 'PD', ?, ?, ?)",
(rid, jobid.strip(), proposal, run,
datetime.now().isoformat(), cmd[3], cmd[4])
)
if "Submitted the following SLURM jobs:" in r:
_, jobids = r.split(":")
jobs = []
for jobid in jobids.split(','):
jobs.append((rid,
jobid.strip(),
proposal,
run,
datetime.now().isoformat(),
cmd[3],
cmd[4])
)
c.executemany(
"INSERT INTO jobs VALUES (?, ?, ?, ?, 'PD', ?, ?, ?)",
jobs)
job_db.commit()
else: # mode == "sim"
......
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