Skip to content
Snippets Groups Projects

[Gotthard-II] Dark Processing.

Merged [Gotthard-II] Dark Processing.
feat/gotthard2_dark into master
Merged Karim Ahmed requested to merge feat/gotthard2_dark into master
Viewing commit 62d7dc3c
Show latest version
2 files
+ 104
− 104
Compare changes
  • Side-by-side
  • Inline
Files
2
notebooks/Gotthard2/Correction_Gotthard2_NBC.ipynb
+ 91
− 98
%% Cell type:markdown id:bed7bd15-21d9-4735-82c1-c27c1a5e3346 tags:
# Gotthard2 Offline Correction #
Author: European XFEL Detector Group, Version: 1.0
Offline Calibration for the Gothard2 Detector
%% Cell type:code id:570322ed-f611-4fd1-b2ec-c12c13d55843 tags:
``` python
in_folder = "/gpfs/exfel/exp/FXE/202221/p003225/raw" # the folder to read data from, required
out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/gotthard2" # the folder to output to, required
run = 50 # run to process, required
sequences = [-1] # sequences to correct, set to [-1] for all, range allowed
sequences_per_node = 1 # number of sequence files per node if notebook executed through xfel-calibrate, set to 0 to not run SLURM parallel
# Parameters used to access raw data.
karabo_id = "FXE_XAD_G2XES" # karabo prefix of Gotthard-II devices
karabo_da = ["GH201"] # data aggregators
receiver_template = "RECEIVER" # receiver template used to read INSTRUMENT keys.
control_template = "CONTROL" # control template used to read CONTROL keys.
instrument_source_template = "{}/DET/{}:daqOutput" # template for source name (filled with karabo_id & receiver_id). e.g. 'SPB_IRDA_JF4M/DET/JNGFR01:daqOutput'
ctrl_source_template = "{}/DET/{}" # template for control source name (filled with karabo_id_control)
karabo_id_control = "" # Control karabo ID. Set to empty string to use the karabo-id
# Parameters for calibration database.
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_interface = "tcp://max-exfl017:8017#8025" # the database interface to use
cal_db_timeout = 180000 # timeout on caldb requests
# Parameters affecting corrected data.
gain_correction = True # do relative gain correction
constants_file = "/gpfs/exfel/data/scratch/ahmedk/dont_remove/gotthard2/constants/calibration_constants_GH2.h5"
# Parameter conditions.
bias_voltage = -1 # Detector bias voltage, set to -1 to use value in raw file.
exposure_time = -1 # Detector exposure time, set to -1 to use value in raw file.
exposure_period = -1 # Detector exposure period, set to -1 to use value in raw file.
operation_mode = -1 # Detector operation mode (1.1/4.5), set to -1 to use value in raw file.
single_photon = -1 # Detector single photon mode (High/Low CDS), set to -1 to use value in raw file.
# Parameters for plotting
skip_plots = False # exit after writing corrected files
pulse_idx_preview = 3 # pulse index to preview. The following even/odd pulse index is used for preview. # TODO: update to pulseId preview.
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:6e9730d8-3908-41d7-abe2-d78e046d5de2 tags:
``` python
import warnings
from functools import partial
import h5py
import pasha as psh
import multiprocessing
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Markdown, display
from extra_data import RunDirectory, H5File
from pathlib import Path
from cal_tools import h5_copy_except
from cal_tools.gotthard2 import gotthard2algs, gotthard2lib
from cal_tools.step_timing import StepTimer
from cal_tools.tools import (
get_constant_from_db_and_time,
get_dir_creation_date,
)
from iCalibrationDB import Conditions, Constants
warnings.filterwarnings('ignore')
prettyPlotting = True
%matplotlib inline
```
%% Cell type:code id:d7c02c48-4429-42ea-a42e-de45366d7fa3 tags:
``` python
in_folder = Path(in_folder)
run_folder = in_folder / f"r{run:04d}"
out_folder = Path(out_folder)
out_folder.mkdir(parents=True, exist_ok=True)
if not karabo_id_control:
karabo_id_control = karabo_id
instrument_src = instrument_source_template.format(
karabo_id, receiver_template)
ctrl_src = ctrl_source_template.format(
karabo_id_control, control_template)
ctrl_src = ctrl_source_template.format(karabo_id_control, control_template)
print(f"Process modules: {karabo_da} for run {run}")
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 not karabo_id_control:
karabo_id_control = karabo_id
# TODO: Remove later
import datetime
creation_time=datetime.datetime.strptime(
"2022-06-28 13:00:00.00",
"%Y-%m-%d %H:%M:%S.%f")
```
%% Cell type:code id:b5eb816e-b5f2-44ce-9907-0273d82341b6 tags:
``` python
# Select only sequence files to process for the selected detector.
if sequences == [-1]:
possible_patterns = list(f"*{mod}*.h5" for mod in karabo_da)
else:
possible_patterns = list(
f"*{mod}-S{s:05d}.h5" for mod in karabo_da for s in sequences)
run_folder = Path(in_folder / f"r{run:04d}")
seq_files = [f for f in run_folder.glob("*.h5") if any(
f.match(p) for p in possible_patterns)]
seq_files = sorted(seq_files)
if not seq_files:
raise IndexError("No sequence files available for the selected sequences.")
print(f"Processing a total of {len(seq_files)} sequence files")
```
%% Cell type:code id:f9a8d1eb-ce6a-4ed0-abf4-4a6029734672 tags:
``` python
step_timer = StepTimer()
```
%% Cell type:code id:892172d8 tags:
``` python
# # Read slow data
# run_dc = RunDirectory(run_folder)
# g2ctrl = gotthard2lib.Gotthard2Ctrl(run_dc=run_dc, ctrl_src=ctrl_src)
# bias_voltage = g2ctrl.get_bias_voltage()
# exposure_time = g2ctrl.get_exposure_time()
# exposure_period = g2ctrl.get_exposure_period()
# operation_mode = g2ctrl.get_operation_mode()
# single_photon = g2ctrl.get_single_photon()
```
%% Cell type:code id:4a8c1b95 tags:
``` python
# condition = Conditions.Dark.gotthard2(
# bias_voltage=bias_voltage,
# exposure_time=exposure_time,
# exposure_period=exposure_period,
# single_photon=single_photon,
# operation_mode=operation_mode,
# )
# stripes = 1280
# def get_constants_for_module(mod: str):
# """Get calibration constants for given module for Gotthard-II."""
# when = {}
# retrieval_function = partial(get_constant_from_db_and_time(
# karabo_id=karabo_id,
# karabo_da=mod,
# cal_db_interface=cal_db_interface,
# creation_time=creation_time,
# timeout=cal_db_timeout,
# print_once=False,
# ))
# constants = {
# "LUT":np.uint32,
# "Offset":np.float32,
# "Noise":np.float32,
# "BadPixelsDark":np.float32
# }
# cons_data[mod] = {}
# for cname, ctype in constants.items():
# cons_data[mod][cname], when[cname] = retrieval_function(
# condition=condition,
# constant=getattr(Constants.gotthard2, cname)(),
# empty_constant=np.zeros((stripes, 2, 3), dtype=ctype)
# )
# if gain_correction:
# for cname in ["BadPixelsFF", "RelativeGain"]:
# cons_data[mod][cname], when[cname] = retrieval_function(
# condition=condition,
# constant=Constants.gotthard2.BadPixelsFF(),
# empty_constant=None
# )
# # combine masks
# if cons_data[mod].get("BadPixelsFF", None) is not None:
# mask |= np.moveaxis(cons_data[mod]["BadPixelsFF"], 0, 1)
# cons_data = {}
# with multiprocessing.Pool() as pool:
# r = pool.map(get_constants_for_module, karabo_da)
# Read slow data
run_dc = RunDirectory(run_folder)
g2ctrl = gotthard2lib.Gotthard2Ctrl(
run_dc=run_dc, ctrl_src=ctrl_src)
bias_voltage = g2ctrl.get_bias_voltage()
exposure_time = g2ctrl.get_exposure_time()
exposure_period = g2ctrl.get_exposure_period()
operation_mode = g2ctrl.get_operation_mode()
single_photon = g2ctrl.get_single_photon()
```
%% Cell type:markdown id:8c852392-bb19-4c40-b2ce-3b787538a92d tags:
### Retrieving calibration constants
%% Cell type:code id:3ba5d81b-e0d6-447b-a7e1-7f40e77b5a61 tags:
%% Cell type:code id:1cdbe818 tags:
``` python
# TODO: Retrieve calibration constants from CALCAT
condition = Conditions.Dark.Gotthard2(
bias_voltage=bias_voltage,
exposure_time=exposure_time,
exposure_period=exposure_period,
# TODO: put back with new data.
# And decide what to do with old data missing these parameters.
# single_photon=single_photon,
# operation_mode=operation_mode,
)
# TODO: Remove temporary constants load.
constants = {}
# load constants temporarily using defined local paths.
with h5py.File(constants_file, 'r') as cfile:
offset_map = cfile["offset_map"][()]
gain_map = cfile["gain_map"][()]
lut = cfile["LUT"][()]
bpix_map = cfile["bpix_ff"][()]
constants[karabo_da[0]] = (lut, offset_map, bpix_map, gain_map)
def get_constants_for_module(mod: str):
"""Get calibration constants for given module for Gotthard-II."""
when = {}
retrieval_function = partial(
get_constant_from_db_and_time,
karabo_id=karabo_id,
karabo_da=mod,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout,
print_once=False,
)
constants = {
"LUT": ((1280, 2, 4096), np.uint16), # TODO: what to do if LUT was not retrieved.
"Offset": ((1280, 2, 3), np.float32),
"Noise": ((1280, 2, 3), np.float32),
"BadPixelsDark": ((1280, 2, 3), np.uint32),
}
cons_data = {}
for cname, cempty in constants.items():
cons_data[cname], when[cname] = retrieval_function(
condition=condition,
constant=getattr(Constants.Gotthard2, cname)(),
empty_constant=np.zeros(cempty[0], dtype=cempty[1])
)
bpix = cons_data["BadPixelsDark"]
if gain_correction:
for cname in ["BadPixelsFF", "RelativeGain"]:
cons_data[cname], when[cname] = retrieval_function(
condition=condition,
constant=getattr(Constants.Gotthard2, cname)(),
empty_constant=None
)
# combine masks
if cons_data.get("BadPixelsFF", None) is not None:
bpix |= cons_data["BadPixelsFF"]
return (
cons_data["LUT"], cons_data["Offset"],
bpix, cons_data["RelativeGain"], mod, when,
)
with multiprocessing.Pool() as pool:
r = pool.map(get_constants_for_module, karabo_da)
```
%% Cell type:code id:8ddbc1ba tags:
``` python
# Print timestamps for the retrieved constants.
# constants = {}
# for offset_map, mask, gain_map, mod, when in r:
# print(f'Constants for module {mod}:')
# for const in when:
# print(f' {const} injected at {when[const]}')
# if gain_map is None:
# print("No gain map found")
# gain_correction = False
# constants[mod] = (offset_map, mask, gain_map)
constants = {}
for lut_map, offset_map, bpix, gain_map, mod, when in r:
print(f'Constants for module {mod}:')
for const in when:
print(f' {const} injected at {when[const]}')
if gain_map is None:
print("No gain map found")
gain_correction = False
constants[mod] = (lut_map, offset_map, bpix, gain_map)
```
%% Cell type:code id:23fcf7f4-351a-4df7-8829-d8497d94fecc tags:
``` python
context = psh.context.ProcessContext(num_workers=multiprocessing.cpu_count())
```
%% Cell type:code id:daecd662-26d2-4cb8-aa70-383a579cf9f9 tags:
``` python
def correct_train(wid, index, d):
g = gain[index]
gotthard2algs.convert_to_10bit(d, lut, data_corr[index, ...])
gotthard2algs.correct_train(
data_corr[index, ...], mask[index, ...], g,
offset_map.astype(np.float32),
gain_map.astype(np.float32),
bpix_map,
offset_map,
gain_map,
bpix,
apply_gain=gain_correction,
)
```
%% Cell type:code id:f88c1aa6-a735-4b72-adce-b30162f5daea tags:
``` python
for mod in karabo_da:
# This is used in case receiver template consists of
# karabo data aggregator index. e.g. detector at DETLAB
instr_mod_src = instrument_src.format(mod[-2:])
data_path = "INSTRUMENT/"+instr_mod_src+"/data"
for raw_file in seq_files:
step_timer.start()
dc = H5File(raw_file)
out_file = out_folder / raw_file.name.replace("RAW", "CORR")
# Select module INSTRUMENT source and deselect empty trains.
dc = dc.select(instr_mod_src, require_all=True)
data = dc[instr_mod_src, "data.adc"].ndarray()
gain = dc[instr_mod_src, "data.gain"].ndarray()
step_timer.done_step("preparing raw data")
dshape = data.shape
step_timer.start()
lut, offset_map, mask, gain_map = constants[mod]
lut, offset_map, bpix, gain_map = constants[mod]
# Allocate shared arrays.
data_corr = context.alloc(shape=dshape, dtype=np.float32)
mask = context.alloc(shape=dshape, dtype=np.uint32)
context.map(correct_train, data)
step_timer.done_step("Correcting one sequence file")
step_timer.start()
# Provided PSI gain map has 0 values. Set inf values to nan.
data_corr[np.isinf(data_corr)] = np.nan
# Create CORR files and add corrected data sources.
# Exclude raw data images (data/adc)
with h5py.File(out_file, 'w') as ofile:
# Copy RAW non-calibrated sources.
with h5py.File(raw_file, 'r') as sfile:
h5_copy_except.h5_copy_except_paths(
sfile, ofile, [f"{data_path}/adc"])
# Create datasets with the available corrected data
ddset = ofile.create_dataset(
f"{data_path}/adc",
data=data_corr,
chunks=((1,) + dshape[1:]), # 1 chunk == 1 image
dtype=np.float32,
)
# Create datasets with the available corrected data
ddset = ofile.create_dataset(
f"{data_path}/mask",
data=mask,
chunks=((1,) + dshape[1:]), # 1 chunk == 1 image
dtype=np.uint32,
compression="gzip",
compression_opts=1,
shuffle=True,
)
step_timer.done_step("Storing data")
```
%% Cell type:code id:94b8e4d2-9f8c-4c23-a509-39238dd8435c tags:
``` python
print(f"Total processing time {step_timer.timespan():.01f} s")
step_timer.print_summary()
```
%% Cell type:code id:0ccc7f7e-2a3f-4ac0-b854-7d505410d2fd tags:
``` python
if skip_plots:
print('Skipping plots')
import sys
sys.exit(0)
```
%% Cell type:code id:ff203f77-3811-46f3-bf7d-226d2dcab13f tags:
``` python
mod_dcs = {}
for mod in karabo_da:
mod_dcs[mod] = {}
with RunDirectory(out_folder) as out_dc:
tid, mod_dcs[mod]["train_corr_data"] = next(
out_dc[instr_mod_src, "data.adc"].trains())
with RunDirectory(run_folder) as in_dc:
train_dict = in_dc.train_from_id(tid)[1][instr_mod_src]
mod_dcs[mod]["train_raw_data"] = train_dict["data.adc"]
mod_dcs[mod]["train_raw_gain"] = train_dict["data.gain"]
```
%% Cell type:code id:1b379438-eb1d-42b2-ac83-eb8cf88c46db tags:
``` python
display(Markdown("### Mean RAW across pulses for one train:"))
display(Markdown(f"Train: {tid}"))
step_timer.start()
for mod in karabo_da:
raw_data = mod_dcs[mod]["train_raw_data"]
fig, ax = plt.subplots(figsize=(20, 10))
im = ax.plot(np.mean(raw_data, axis=0))
ax.set_title(f'Module {mod}')
ax.set_xlabel("Strip #", size=15)
ax.set_ylabel("12-bit ADC output", size=15)
step_timer.done_step("Plotting raw data")
```
%% Cell type:code id:512b48a6-9e2d-4fcc-b075-3698297177bf tags:
``` python
display(Markdown("### Mean CORRECTED across pulses for one train:"))
display(Markdown(f"Train: {tid}"))
step_timer.start()
for mod in karabo_da:
corr_data = mod_dcs[mod]["train_corr_data"]
fig, ax = plt.subplots(figsize=(20, 10))
im = ax.plot(np.mean(corr_data, axis=0))
ax.set_title(f'Module {mod}')
ax.set_xlabel("Strip #", size=20)
ax.set_ylabel("10-bit KeV. output", size=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
step_timer.done_step("Plotting corrected data")
```
%% Cell type:code id:cd8f5e08-fcee-4bff-ba63-6452b3d892a2 tags:
``` python
# Validate given "pulse_idx_preview"
if pulse_idx_preview + 1 > data.shape[1]:
print(f"WARNING: selected pulse_idx_preview {pulse_idx_preview} is not available in data."
" Previewing 1st pulse.")
pulse_idx_preview = 1
if data.shape[1] == 1:
odd_pulse = 1
even_pulse = None
else:
odd_pulse = pulse_idx_preview if pulse_idx_preview % 2 else pulse_idx_preview + 1
even_pulse = pulse_idx_preview if not (pulse_idx_preview % 2) else pulse_idx_preview + 1
if pulse_idx_preview + 1 > data.shape[1]:
pulse_idx_preview = 1
if data.shape[1] > 1:
pulse_idx_preview = 2
```
%% Cell type:code id:e5f0d4d8-e32c-4f2c-8469-4ebbfd3f644c tags:
``` python
display(Markdown("### RAW even/odd pulses for one train:"))
display(Markdown(f"RAW train: {tid}"))
for mod in karabo_da:
raw_data = mod_dcs[mod]["train_raw_data"]
fig, ax = plt.subplots(figsize=(20, 10))
ax.plot(raw_data[odd_pulse], label=f"Odd Pulse {odd_pulse}")
if even_pulse:
ax.plot(raw_data[even_pulse], label=f"Even Pulse {even_pulse}")
ax.set_title(f'Module {mod}')
ax.set_xlabel("Strip #", size=20)
ax.set_ylabel("12-bit ADC output", size=20)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
ax.legend()
step_timer.done_step("Plotting RAW odd/even pulses.")
```
%% Cell type:markdown id:f1ddac79-5590-4bdb-9cfa-147d73f5dd5f tags:
### Single RAW Train Preview
A single image from the first RAW train
%% Cell type:code id:96900ba8-96a3-44b0-b489-d948f19298dc tags:
``` python
display(Markdown("### CORRECTED even/odd pulses for one train:"))
display(Markdown(f"CORRECTED train: {tid}"))
for mod in karabo_da:
corr_data = mod_dcs[mod]["train_corr_data"]
fig, ax = plt.subplots(figsize=(20, 10))
ax.plot(corr_data[odd_pulse], label=f"Odd Pulse {odd_pulse}")
if even_pulse:
ax.plot(corr_data[even_pulse], label=f"Even Pulse {even_pulse}")
ax.set_title(f'Module {mod}')
ax.set_xlabel("Strip #", size=15)
ax.set_ylabel("10-bit KeV. output", size=15)
ax.legend()
step_timer.done_step("Plotting CORRECTED odd/even pulses.")
```

Legal Notice | Contacts