Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • calibration/pycalibration
1 result
Show changes
Commits on Source (19)
%% Cell type:markdown id: tags:
# Characterization of dark and flat field for Dynamic Flat Field correction
Author: Egor Sobolev
Computation of dark offsets and flat-field principal components
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/SPB/202430/p900425/raw" # input folder, required
out_folder = '/gpfs/exfel/data/scratch/esobolev/test/shimadzu' # output folder, required
metadata_folder = "" # Directory containing calibration_metadata.yml when run by xfel-calibrate
run_high = 1 # run number in which dark data was recorded, required
run_low = 2 # run number in which flat-field data was recorded, required
operation_mode = "PCA_DynamicFF" # Detector operation mode, optional (defaults to "PCA_DynamicFF")
# Data files parameters.
karabo_da = ['-1'] # data aggregators
karabo_id = "SPB_MIC_HPVX2" # karabo prefix of Shimadzu HPV-X2 devices
# Database access parameters.
cal_db_interface = "tcp://max-exfl-cal001:8021" # Unused, calibration DB interface to use
cal_db_timeout = 30000 # Unused, calibration DB timeout
db_output = False # if True, the notebook sends dark constants to the calibration database
local_output = True # if True, the notebook saves dark constants locally
creation_time = "" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HH:MM:SS.00 e.g. 2019-07-04 11:02:41.00
# Calibration constants parameters
frame_range = [4] # range list [start, end, step] of frame indices within a train to use for characterization.
n_components = 50 # Number of principal components of flat-field to compute (default: 50)
```
%% Cell type:code id: tags:
``` python
import datetime
import os
import warnings
from logging import warning
from shutil import copyfile
from tempfile import NamedTemporaryFile
warnings.filterwarnings('ignore')
import time
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, Markdown
from extra_data import RunDirectory
%matplotlib inline
from cal_tools.step_timing import StepTimer
from cal_tools.tools import (
get_dir_creation_date,
run_prop_seq_from_path,
save_dict_to_hdf5
calcat_creation_time,
)
from cal_tools.restful_config import calibration_client, extra_calibration_client
from cal_tools.shimadzu import ShimadzuHPVX2
from cal_tools.constants import write_ccv, inject_ccv
from cal_tools.constants import CCVAlreadyInjectedError, write_ccv, inject_ccv
import dynflatfield as dffc
from dynflatfield.draw import plot_images, plot_camera_image
```
%% Cell type:code id: tags:
``` python
extra_calibration_client() # Configure CalibrationData.
cc = calibration_client()
pdus = cc.get_all_phy_det_units_from_detector(
{"detector_identifier": karabo_id}) # TODO: Use creation_time for snapshot_at
pdus = cc.get_all_phy_det_units_from_detector({
"detector_identifier": karabo_id,
"pdu_snapshot_at": calcat_creation_time(
in_folder, min(run_low, run_high), creation_time),
})
if not pdus["success"]:
raise ValueError("Failed to retrieve PDUs")
detector_info = pdus['data'][0]['detector']
detector = ShimadzuHPVX2(detector_info["source_name_pattern"])
print(f"Instrument {detector.instrument}")
print(f"Detector in use is {karabo_id}")
modules = {}
for pdu_no, pdu in enumerate(pdus["data"]):
db_module = pdu["physical_name"]
module = pdu["module_number"]
da = pdu["karabo_da"]
if karabo_da[0] != "-1" and da not in karabo_da:
continue
instrument_source_name = detector.instrument_source(module)
print('-', da, db_module, module, instrument_source_name)
modules[da] = dict(
db_module=db_module,
module=module,
raw_source_name=instrument_source_name,
pdu_no=pdu_no,
)
constants = {}
# make a frame slice
frame_range += [None] * (3 - len(frame_range))
frame_slice = slice(*frame_range)
step_timer = StepTimer()
```
%% Cell type:markdown id: tags:
# Offset map
%% Cell type:code id: tags:
``` python
dark_run = run_high
dark_creation_time = get_dir_creation_date(in_folder, dark_run)
dark_creation_time = calcat_creation_time(in_folder, dark_run, creation_time)
print(f"Using {dark_creation_time} as creation time of Offset constant.")
for da, meta in modules.items():
source_name = detector.instrument_source(meta["module"])
image_key = detector.image_key
display(Markdown(f"## {source_name}"))
# read
step_timer.start()
file_da, _, _ = da.partition('/')
dark_dc = RunDirectory(f"{in_folder}/r{dark_run:04d}",
include=f"RAW-R{dark_run:04d}-{file_da}-S*.h5")
if source_name not in dark_dc.all_sources:
raise ValueError(f"Could not find source {source_name} for module {da} in dark data")
dark_dc = dark_dc.select([(source_name, image_key)])
conditions = detector.conditions(dark_dc, meta["module"])
key_data = dark_dc[source_name, image_key]
images_dark = key_data.ndarray()
images_dark = key_data.ndarray(roi=frame_slice)
ntrain, npulse, ny, nx = images_dark.shape
print(f"N image: {ntrain * npulse} (ntrain: {ntrain}, npulse: {npulse})")
print(f"N image: {ntrain * npulse} (ntrain: {ntrain}, npulse: {npulse}/{key_data.shape[1]})")
print(f"Image size: {ny} x {nx} px")
step_timer.done_step("Read dark images")
# process
step_timer.start()
dark = dffc.process_dark(images_dark) # Amounts to a per-pixel mean right now.
# put results in the dict
module_constants = constants.setdefault(meta["db_module"], {})
module_constants["Offset"] = dict(
conditions=conditions, data=dark, pdu_no=meta["pdu_no"],
creation_time=dark_creation_time, dims=['ss', 'fs']
)
step_timer.done_step("Process dark images")
display()
# draw plots
step_timer.start()
plot_camera_image(dark)
plt.show()
step_timer.done_step("Draw offsets")
```
%% Cell type:markdown id: tags:
# Flat-field PCA decomposition
%% Cell type:code id: tags:
``` python
flat_run = run_low
flat_creation_time = get_dir_creation_date(in_folder, flat_run)
flat_creation_time = calcat_creation_time(in_folder, flat_run, creation_time)
print(f"Using {flat_creation_time} as creation time of DynamicFF constant.")
for da, meta in modules.items():
source_name = detector.instrument_source(meta["module"])
image_key = detector.image_key
display(Markdown(f"## {source_name}"))
# read
step_timer.start()
file_da, _, _ = da.partition('/')
flat_dc = RunDirectory(f"{in_folder}/r{flat_run:04d}",
include=f"RAW-R{flat_run:04d}-{file_da}-S*.h5")
if source_name not in flat_dc.all_sources:
raise ValueError(f"Could not find source {source_name} for module {da} in flatfield data")
flat_dc = flat_dc.select([(source_name, image_key)])
conditions = detector.conditions(flat_dc, meta["module"])
dark = constants[meta["db_module"]]["Offset"]["data"]
dark_conditions = constants[meta["db_module"]]["Offset"]["conditions"]
if conditions != dark_conditions:
raise ValueError(f"The conditions for flat-field run {conditions}) do not match "
f"the dark run conditions ({dark_conditions}). Skip flat-field characterization.")
key_data = flat_dc[source_name][image_key]
images_flat = key_data.ndarray()
images_flat = key_data.ndarray(roi=frame_slice)
ntrain, npulse, ny, nx = images_flat.shape
print(f"N image: {ntrain * npulse} (ntrain: {ntrain}, npulse: {npulse})")
print(f"N image: {ntrain * npulse} (ntrain: {ntrain}, npulse: {npulse}/{key_data.shape[1]})")
print(f"Image size: {ny} x {nx} px")
step_timer.done_step("Read flat-field images")
# process
step_timer.start()
flat, components, explained_variance_ratio = dffc.process_flat(
images_flat, dark, n_components)
flat_data = np.concatenate([flat[None, ...], components])
# put results in the dict
conditions = detector.conditions(flat_dc, meta["module"])
module_constants = constants.setdefault(meta["db_module"], {})
module_constants["DynamicFF"] = dict(
conditions=conditions, data=flat_data, pdu_no=meta["pdu_no"],
creation_time=flat_creation_time, dims=['component', 'ss', 'fs']
)
step_timer.done_step("Process flat-field images")
# draw plots
step_timer.start()
display(Markdown("### Average flat-field"))
plot_camera_image(flat)
plt.show()
display(Markdown("### Explained variance ratio"))
fig, ax = plt.subplots(1, 1, figsize=(10,4), tight_layout=True)
ax.semilogy(explained_variance_ratio, 'o')
ax.set_xticks(np.arange(len(explained_variance_ratio)))
ax.set_xlabel("Component no.")
ax.set_ylabel("Variance fraction")
plt.show()
display(Markdown("### The first principal components (up to 20)"))
plot_images(components[:20], figsize=(13, 8))
plt.show()
step_timer.done_step("Draw flat-field")
```
%% Cell type:markdown id: tags:
## Calibration constants
%% Cell type:code id: tags:
``` python
step_timer.start()
_, proposal, _ = run_prop_seq_from_path(in_folder)
# Output Folder Creation:
if local_output:
os.makedirs(out_folder, exist_ok=True)
for db_module, module_constants in constants.items():
for constant_name, constant in module_constants.items():
conditions = constant["conditions"]
pdu = pdus["data"][constant["pdu_no"]]
with NamedTemporaryFile() as tempf:
ccv_root = write_ccv(
tempf.name,
pdu['physical_name'], pdu['uuid'], pdu['detector_type']['name'],
constant_name, conditions, constant['creation_time'],
proposal, [dark_run, flat_run],
constant["data"], constant['dims'])
if db_output:
inject_ccv(tempf.name, ccv_root, metadata_folder)
try:
inject_ccv(tempf.name, ccv_root, metadata_folder)
except CCVAlreadyInjectedError:
pass
if local_output:
ofile = f"{out_folder}/const_{constant_name}_{db_module}.h5"
if os.path.isfile(ofile):
print(f'File {ofile} already exists and will be overwritten')
copyfile(tempf.name, ofile)
```
%% Cell type:code id: tags:
``` python
print(f"Total processing time {step_timer.timespan():.01f} s")
step_timer.print_summary()
```
......
%% Cell type:markdown id: tags:
# Dynamic Flat-field Offline Correction
Author: Egor Sobolev
Offline dynamic flat-field correction
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/SPB/202430/p900425/raw" # input folder, required
out_folder ="/gpfs/exfel/exp/SPB/202430/p900425/scratch/proc/r0003" # output folder, required
metadata_folder = "" # Directory containing calibration_metadata.yml when run by xfel-calibrate
run = 3 # which run to read data from, required
# Data files parameters.
karabo_da = ['-1'] # data aggregators
karabo_id = "SPB_MIC_HPVX2" # karabo prefix of Shimadzu HPV-X2 devices
# Database access parameters.
cal_db_interface = "tcp://max-exfl-cal001:8021" # Unused, calibration DB interface to use
cal_db_timeout = 30000 # Unused, calibration DB timeout
creation_time = "" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HH:MM:SS.00 e.g. 2019-07-04 11:02:41.00
# Correction parameters
n_components = 20 # number of principal components of flat-field to use in correction
downsample_factors = [1, 1] # list of downsample factors for each image dimention (y, x)
num_proc = 32 # number of processes running correction in parallel
```
%% Cell type:code id: tags:
``` python
import os
import sys
import h5py
import warnings
from logging import warning
warnings.filterwarnings('ignore')
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, Markdown
from datetime import datetime
from extra_data import RunDirectory, by_id
%matplotlib inline
from cal_tools.step_timing import StepTimer
from cal_tools.files import sequence_trains, DataFile
from cal_tools.tools import get_dir_creation_date
from cal_tools.tools import calcat_creation_time
from cal_tools.restful_config import calibration_client, extra_calibration_client
from cal_tools.calcat_interface2 import CalibrationData
from cal_tools.shimadzu import ShimadzuHPVX2
from dynflatfield import (
DynamicFlatFieldCorrectionCython as DynamicFlatFieldCorrection,
FlatFieldCorrectionFileProcessor
FileDynamicFlatFieldProcessor
)
from dynflatfield.draw import plot_images, plot_camera_image
```
%% Cell type:code id: tags:
``` python
creation_time = get_dir_creation_date(in_folder, run)
creation_time = calcat_creation_time(in_folder, run, creation_time)
print(f"Creation time is {creation_time}")
extra_calibration_client() # Configure CalibrationData API.
cc = calibration_client()
pdus = cc.get_all_phy_det_units_from_detector(
{"detector_identifier": karabo_id}) # TODO: Use creation_time for snapshot_at
pdus = cc.get_all_phy_det_units_from_detector({
"detector_identifier": karabo_id,
"pdu_snapshot_at": creation_time,
})
if not pdus["success"]:
raise ValueError("Failed to retrieve PDUs")
detector_info = pdus['data'][0]['detector']
detector = ShimadzuHPVX2(detector_info["source_name_pattern"])
index_group = detector.image_index_group
image_key = detector.image_key
print(f"Instrument {detector.instrument}")
print(f"Detector in use is {karabo_id}")
modules = {}
for pdu in pdus["data"]:
db_module = pdu["physical_name"]
module = pdu["module_number"]
da = pdu["karabo_da"]
if karabo_da[0] != "-1" and da not in karabo_da:
continue
instrument_source_name = detector.instrument_source(module)
corrected_source_name = detector.corrected_source(module)
print('-', da, db_module, module, instrument_source_name)
modules[da] = dict(
db_module=db_module,
module=module,
raw_source_name=instrument_source_name,
corrected_source_name=corrected_source_name,
)
step_timer = StepTimer()
```
%% Cell type:markdown id: tags:
# Calibration constants
%% Cell type:code id: tags:
``` python
step_timer.start()
dc = RunDirectory(f"{in_folder}/r{run:04d}")
conditions = detector.conditions(dc)
caldata = CalibrationData.from_condition(
conditions, 'SPB_MIC_HPVX2', event_at=creation_time)
aggregators = {}
corrections = {}
for da in modules:
file_da, _, _ = da.partition('/')
aggregators.setdefault(file_da, []).append(da)
try:
dark = caldata["Offset", da].ndarray()
flat = caldata["DynamicFF", da].ndarray()
components = flat[1:][:n_components]
flat = flat[0]
dffc = DynamicFlatFieldCorrection.from_constants(
dark, flat, components, downsample_factors)
corrections[da] = dffc
file_da, _, _ = da.partition('/')
aggregators.setdefault(file_da, []).append(da)
except (KeyError, FileNotFoundError):
warning(f"Constants are not found for module {da}. "
"The module will not calibrated")
# missed constants are reported later
pass
step_timer.done_step("Load calibration constants")
```
%% Cell type:markdown id: tags:
# Correction
%% Cell type:code id: tags:
``` python
# Output Folder Creation:
os.makedirs(out_folder, exist_ok=True)
report = []
missed_constants = 0
report = {}
for file_da, file_modules in aggregators.items():
dc = RunDirectory(f"{in_folder}/r{run:04d}", f"RAW-R{run:04d}-{file_da}-S*.h5")
# build train IDs
train_ids = set()
process_modules = []
for da in file_modules:
instrument_source = modules[da]["raw_source_name"]
if instrument_source in dc.all_sources:
keydata = dc[instrument_source][image_key].drop_empty_trains()
train_ids.update(keydata.train_ids)
process_modules.append(da)
else:
if instrument_source not in dc.all_sources:
print(f"Source {instrument_source} for module {da} is missed")
continue
if da not in corrections:
missed_constants += 1
warning(f"Constants are not found for module {da}. "
"The module will not calibrated")
continue
keydata = dc[instrument_source][image_key].drop_empty_trains()
train_ids.update(keydata.train_ids)
process_modules.append(da)
if not process_modules:
continue
train_ids = np.array(sorted(train_ids))
ts = dc.select_trains(by_id[train_ids]).train_timestamps().astype(np.uint64)
# correct and write sequence files
seq_report = {}
for seq_id, train_mask in sequence_trains(train_ids, 200):
step_timer.start()
print('* sequence', seq_id)
seq_train_ids = train_ids[train_mask]
seq_timestamps = ts[train_mask]
dc_seq = dc.select_trains(by_id[seq_train_ids])
ntrains = len(seq_train_ids)
# create output file
channels = [f"{modules[da]['corrected_source_name']}/{index_group}"
for da in process_modules]
f = DataFile.from_details(out_folder, file_da, run, seq_id)
f.create_metadata(like=dc, instrument_channels=channels)
f.create_index(seq_train_ids, timestamps=seq_timestamps)
# create file structure
seq_report = {}
file_datasets = {}
for da in process_modules:
instrument_source = modules[da]["raw_source_name"]
keydata = dc_seq[instrument_source][image_key].drop_empty_trains()
count = keydata.data_counts(labelled=False)
i = np.flatnonzero(count)
raw_images = keydata.select_trains(np.s_[i]).ndarray()
# not pulse resolved
shape = keydata.shape
count = np.in1d(seq_train_ids, keydata.train_ids).astype(int)
corrected_source = modules[da]["corrected_source_name"]
src = f.create_instrument_source(corrected_source)
src.create_index(**{index_group: count})
# create key for images
ds_data = src.create_key(image_key, shape=shape, dtype=np.float32)
module_datasets = {image_key: ds_data}
# create keys for image parameters
for key in detector.copy_keys:
keydata = dc_seq[instrument_source][key].drop_empty_trains()
module_datasets[key] = (keydata, src.create_key(
key, shape=keydata.shape, dtype=keydata.dtype))
file_datasets[da] = module_datasets
step_timer.done_step("Create output file")
# correct and write data to file
for da in process_modules:
step_timer.start()
dc_seq = dc.select_trains(by_id[seq_train_ids])
dffc = corrections[da]
instrument_source = modules[da]["raw_source_name"]
proc = FlatFieldCorrectionFileProcessor(dffc, num_proc, instrument_source, image_key)
proc = FileDynamicFlatFieldProcessor(dffc, num_proc, instrument_source, image_key)
proc.start_workers()
proc.run(dc_seq)
proc.join_workers()
# not pulse resolved
corrected_images = np.stack(proc.rdr.results, 0)
file_datasets[da][image_key][:] = corrected_images
# copy image parameters
for key in detector.copy_keys:
keydata, ds = file_datasets[da][key]
ds[:] = keydata.ndarray()
seq_report[da] = (raw_images[0, 0], corrected_images[:20, 0])
rep_cix = np.argmax(np.mean(raw_images[:20], axis=(-1, -2)), axis=1)
rep_rix = rep_cix[0]
da_report = seq_report.setdefault(da, [])
da_report.append((raw_images[0, rep_rix],
corrected_images[range(len(rep_cix)), rep_cix]))
step_timer.done_step("Correct flat-field")
f.close()
report.append(seq_report)
report.update(seq_report)
if not report:
if missed_constants:
raise ValueError("Calibration constants are not found for any module")
else:
warning("No data to correct")
sys.exit(0)
```
%% Cell type:code id: tags:
``` python
step_timer.start()
if report:
for da, (raw_image, corrected_images) in report[0].items():
for da, da_report in report.items():
if len(da_report) > 0:
raw_images, corrected_images = zip(*da_report)
raw_images = np.stack(raw_images)
corrected_images = np.concatenate(corrected_images, axis=0)
else:
raw_images, corrected_images = da_report
raw_images = raw_images[None, ...]
source = modules[da]["raw_source_name"]
display(Markdown(f"## {source}"))
display(Markdown("### The first raw image"))
plot_camera_image(raw_images[0, 0])
display(Markdown("### The brightest raw image from the first train"))
plot_camera_image(raw_images[0])
plt.show()
display(Markdown("### The first corrected image"))
display(Markdown("### The brightest corrected image from the first train"))
plot_camera_image(corrected_images[0])
plt.show()
min_images = 5
if corrected_images.shape[0] < min_images:
# Update corrected_images to avoid less axes
# array dimension than expected in plot_images.
corrected_images = np.concatenate(
[
corrected_images,
np.full(
(min_images - corrected_images.shape[0], *corrected_images.shape[1:]),
fill_value=np.nan)])
display(Markdown("### The first corrected images in the trains (up to 20)"))
display(Markdown("### The brightest corrected images by one from a train (up to 20 first trains)"))
plot_images(corrected_images, figsize=(13, 8))
plt.show()
step_timer.done_step("Draw images")
```
%% Cell type:code id: tags:
``` python
print(f"Total processing time {step_timer.timespan():.01f} s")
step_timer.print_summary()
```
......
%% Cell type:markdown id:49b6577f-96a5-4dd2-bdd9-da661b2c4619 tags:
# Gotthard2 Dark Image Characterization #
Author: European XFEL Detector Group, Version: 1.0
The following is a processing for the dark constants (`Offset`, `Noise`, and `BadPixelsDark`) maps using dark images taken with Gotthard2 detector (GH2 50um or 25um).
All constants are evaluated per strip, per pulse, and per memory cell. The maps are calculated for each gain stage that is acquired in 3 separate runs.
The three maps are of shape (stripes, cells, gains): (1280, 2, 3). They can be injected to the database (`db_output`) and/or stored locally (`local_output`).
%% Cell type:code id:818e24e8 tags:
``` python
in_folder = "/gpfs/exfel/exp/FXE/202231/p900298/raw" # the folder to read data from, required
out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/gotthard2/darks" # the folder to output to, required
metadata_folder = '' # Directory containing calibration_metadata.yml when run by xfel-calibrate
run_high = 7 # run number for G0 dark run, required
run_med = 8 # run number for G1 dark run, required
run_low = 9 # run number for G2 dark run, required
# Parameters used to access raw data.
input_source_template = "{karabo_id}/DET/RECEIVER{input_source_affixes}:daqOutput" # data source template used to read INSTRUMENT keys.
ctrl_source_template = "{karabo_id_control}/DET/{control_template}" # template for control source name (filled with karabo_id_control)
karabo_id = "FXE_XAD_G2XES" # karabo prefix of Gotthard-II devices
karabo_da = [""] # data aggregators
input_source_affixes = [""] # The affix to format into the input source template to be able to load the correct module's data.
control_template = "CONTROL" # control template used to read CONTROL keys.
karabo_id_control = "" # Control karabo ID. Set to empty string to use the karabo-id
# Parameters for the calibration database.
cal_db_interface = "tcp://max-exfl-cal001:8020" # calibration DB interface to use
cal_db_timeout = 300000 # timeout on caldb requests
creation_time = "" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC e.g. "2022-06-28 13:00:00"
db_output = False # Output constants to the calibration database
local_output = True # Output constants locally
# Conditions used for injected calibration constants.
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.
acquisition_rate = -1. # Detector acquisition rate (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.
# This is used for the summary notebook
sensor_type = "" # Specifies the GH2 sensor type: '50um', '25um', or "" to automatically decide from RUN data.
# Parameters used for calculating calibration constants
bpix_noise_nitrs = 5 # number of maximum iterations to actively try to reach correct number of badpixels of NOISE_OUT_OF_THRESHOLD type. Leave -1 to keep iterating until finding maximum number of badpixels.
# Parameters used during selecting raw data trains.
min_trains = 1 # Minimum number of trains that should be available to process dark constants. Default 1.
max_trains = 1000 # Maximum number of trains to use for processing dark constants. Set to 0 to use all available trains.
badpixel_threshold_sigma = 5. # bad pixels defined by values outside n times this std from median
# Don't delete! myMDC sends this by default.
operation_mode = '' # Detector dark run acquiring operation mode, optional
```
%% Cell type:code id:8085f9aa tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import pasha as psh
from IPython.display import Markdown, display
from extra_data import RunDirectory
from pathlib import Path
import yaml
from cal_tools.calcat_interface import CalCatApi
from cal_tools.enums import BadPixels
from cal_tools.gotthard2 import gotthard2algs, gotthard2lib
from cal_tools.step_timing import StepTimer
from cal_tools.restful_config import calibration_client
from cal_tools.tools import (
calcat_creation_time,
get_constant_from_db_and_time,
get_report,
raw_data_location_string,
save_const_to_h5,
send_to_db,
)
from iCalibrationDB import Conditions, Constants
%matplotlib inline
```
%% Cell type:code id:bc4a9824 tags:
``` python
# Input parameters validation
gotthard2lib.validate_sensor_type(sensor_type)
```
%% Cell type:code id:18fe4379 tags:
``` python
run_nums = [run_high, run_med, run_low]
in_folder = Path(in_folder)
out_folder = Path(out_folder)
out_folder.mkdir(parents=True, exist_ok=True)
if not karabo_id_control:
karabo_id_control = karabo_id
ctrl_src = ctrl_source_template.format(
karabo_id_control=karabo_id_control,
control_template=control_template,
)
run_nums = gotthard2lib.sort_dark_runs_by_gain(
raw_folder=in_folder,
runs=run_nums,
ctrl_src=ctrl_src,
)
# Read report path to associate it later with injected constants.
report = get_report(metadata_folder)
# Run's creation time:
creation_time = calcat_creation_time(in_folder, run_nums[0], creation_time)
print(f"Creation time: {creation_time}")
```
%% Cell type:code id:c176a86f tags:
``` python
run_dc = RunDirectory(in_folder / f"r{run_nums[0]:04d}")
proposal = list(filter(None, str(in_folder).strip('/').split('/')))[-2]
file_loc = raw_data_location_string(proposal, run_nums)
```
%% Cell type:code id:108be688 tags:
``` python
step_timer = StepTimer()
```
%% Cell type:code id:ff9149fc tags:
``` python
# Read parameter conditions and validate the values for the three runs.
step_timer.start()
run_dcs_dict = dict()
conditions = {
"bias_voltage": set(),
"exposure_time": set(),
"exposure_period": set(),
"acquisition_rate": set(),
"single_photon": set(),
}
for gain, run in enumerate(run_nums):
run_dc = RunDirectory(in_folder / f"r{run:04d}/")
run_dcs_dict[run] = [gain, run_dc]
# Read slow data.
g2ctrl = gotthard2lib.Gotthard2Ctrl(run_dc=run_dc, ctrl_src=ctrl_src)
conditions["bias_voltage"].add(
g2ctrl.get_bias_voltage() if bias_voltage == -1 else bias_voltage
)
conditions["exposure_time"].add(
g2ctrl.get_exposure_time() if exposure_time == -1 else exposure_time
)
conditions["exposure_period"].add(
g2ctrl.get_exposure_period() if exposure_period == -1 else exposure_period
)
conditions["single_photon"].add(
g2ctrl.get_single_photon() if single_photon == -1 else single_photon
)
conditions["acquisition_rate"].add(
g2ctrl.get_acquisition_rate() if acquisition_rate == -1 else acquisition_rate
)
for c, v in conditions.items():
assert len(v) == 1, f"{c} value is not the same for the three runs."
bias_voltage = conditions["bias_voltage"].pop()
print("Bias voltage: ", bias_voltage)
exposure_time = conditions["exposure_time"].pop()
print("Exposure time: ", exposure_time)
exposure_period = conditions["exposure_period"].pop()
print("Exposure period: ", exposure_period)
single_photon = conditions["single_photon"].pop()
print("Single photon: ", single_photon)
acquisition_rate = conditions["acquisition_rate"].pop()
print("Acquisition rate: ", acquisition_rate)
gh2_detector = g2ctrl.get_det_type()
print(f"Processing {gh2_detector} Gotthard2.")
if not sensor_type:
sensor_type = g2ctrl.get_sensor_type()
print(f"Processing {sensor_type} Gotthard2.")
```
%% Cell type:code id:f64bc150-cfcd-4f98-83f9-a982fdacedd7 tags:
``` python
calcat = CalCatApi(client=calibration_client())
detector_id = calcat.detector(karabo_id)['id']
pdus_by_da = calcat.physical_detector_units(detector_id, pdu_snapshot_at=creation_time)
# Module index is stored in the da_to_pdu dict to be used
# later in selecting the correct data source.
# TODO: This can be improved later by getting the module index and even data source
# when we start storing in the CALCAT.
da_to_pdu = dict()
for i, (da, p) in enumerate(sorted(pdus_by_da.items())):
da_to_pdu[da] = p['physical_name']
calcat_das = list(da_to_pdu.keys())
if karabo_da != [""]:
# Filter DA connected to detector in CALCAT
karabo_da = [da for da in karabo_da if da in calcat_das]
else:
karabo_da = calcat_das
print(f"Processing {karabo_da}")
da_to_source = gotthard2lib.map_da_to_source(
run_dc,
calcat_das,
input_source_template,
karabo_id,
input_source_affixes,
)
```
%% Cell type:code id:ac9c5dc3-bc66-4e7e-b6a1-360259be535c tags:
``` python
def specify_trains_to_process(
img_key_data: "extra_data.KeyData",
run_num: int,
src: str,
):
"""Specify total number of trains to process.
Based on given min_trains and max_trains, if given.
Print number of trains to process and number of empty trains.
Raise ValueError if specified trains are less than min_trains.
"""
# Specifies total number of trains to process.
n_trains = img_key_data.shape[0]
all_trains = len(img_key_data.train_ids)
print(
f"{src} for run {run} has {all_trains - n_trains} "
f"trains with empty frames out of {all_trains} trains"
)
if n_trains < min_trains:
raise ValueError(
f"Less than {min_trains} trains are available in RAW data."
" Not enough data to process darks."
)
if max_trains > 0:
n_trains = min(n_trains, max_trains)
print(f"Processing {n_trains} trains.")
return n_trains
```
%% Cell type:code id:3c59c11d tags:
``` python
# set the operating condition
condition = Conditions.Dark.Gotthard2(
bias_voltage=bias_voltage,
exposure_time=exposure_time,
exposure_period=exposure_period,
acquisition_rate=acquisition_rate,
single_photon=single_photon,
)
```
%% Cell type:code id:e2eb2fc0-df9c-4887-9691-f81474f8c131 tags:
``` python
def convert_train(wid, index, tid, d):
"""Convert a Gotthard2 train from 12bit to 10bit."""
gotthard2algs.convert_to_10bit(
d[src]["data.adc"], lut, data_10bit[index, ...]
)
```
%% Cell type:code id:4e8ffeae tags:
``` python
# Calculate noise and offset per pixel and global average, std and median
noise_map = dict()
offset_map = dict()
badpixels_map = dict()
context = psh.ProcessContext(num_workers=25)
empty_lut = (np.arange(2 ** 12).astype(np.float64) * 2 ** 10 / 2 ** 12).astype(
np.uint16
)
empty_lut = np.stack(1280 * [np.stack([empty_lut] * 2)], axis=0)
for mod in karabo_da:
# For 50um idx always 0. For 25um pick the right data source
# from list of 2 sources.
src = da_to_source[mod]
# Retrieve LUT constant
lut, time = get_constant_from_db_and_time(
constant=Constants.Gotthard2.LUT(),
condition=condition,
empty_constant=empty_lut,
karabo_id=karabo_id,
karabo_da=mod,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout,
print_once=False,
)
print(f"Retrieved LUT constant with creation-time {time}")
cshape = (1280, 2, 3)
offset_map[mod] = context.alloc(shape=cshape, dtype=np.float32)
noise_map[mod] = context.alloc(like=offset_map[mod])
badpixels_map[mod] = context.alloc(like=offset_map[mod], dtype=np.uint32)
for run_num, [gain, run_dc] in run_dcs_dict.items():
step_timer.start()
n_trains = specify_trains_to_process(run_dc[src, "data.adc"], run_num, src)
# Select requested number of trains to process.
dc = run_dc.select(src, require_all=True).select_trains(
np.s_[:n_trains]
)
step_timer.done_step("preparing raw data")
step_timer.start()
# Convert 12bit data to 10bit
data_10bit = context.alloc(
shape=dc[src, "data.adc"].shape, dtype=np.float32
)
context.map(convert_train, dc)
step_timer.done_step("convert to 10bit")
step_timer.start()
# The first ~20 frames of each train are excluded.
# The electronics needs some time to reach stable operational conditions
# at the beginning of each acquisition cycle,
# hence the first ~20 images have lower quality and should not be used.
# The rough number of 20 is correct at 4.5 MHz acquisition rate,
# 5 should be enough at 1.1 MHz. Sticking with 20 excluded frames, as there isn't
# much of expected difference.
# Split even and odd data to calculate the two storage cell constants.
# Detector always operates in burst mode.
even_data = data_10bit[:, 20::2, :]
odd_data = data_10bit[:, 21::2, :]
def offset_noise_cell(wid, index, d):
offset_map[mod][:, index, gain] = np.mean(d, axis=(0, 1))
noise_map[mod][:, index, gain] = np.std(d, axis=(0, 1))
# Calculate Offset and Noise Maps.
context.map(offset_noise_cell, (even_data, odd_data))
# Split even and odd gain data.
data_gain = dc[src, "data.gain"].ndarray()
even_gain = data_gain[:, 20::2, :]
odd_gain = data_gain[:, 21::2, :]
raw_g = 3 if gain == 2 else gain
def badpixels_cell(wid, index, g):
"""Check if there are wrong bad gain values.
Indicate pixels with wrong gain value across all trains for each cell."""
badpixels_map[mod][
np.mean(g, axis=(0, 1)) != raw_g, index, :
] |= BadPixels.WRONG_GAIN_VALUE.value
context.map(badpixels_cell, (even_gain, odd_gain))
step_timer.done_step("Processing darks")
```
%% Cell type:code id:8e410ca2 tags:
``` python
print(f"Total processing time {step_timer.timespan():.01f} s")
step_timer.print_summary()
```
%% Cell type:code id:3fc17e05-17ab-4ac4-9e79-c95399278bb9 tags:
``` python
def print_bp_entry(bp):
print(f"{bp.name:<30s} {bp.value:032b} -> {int(bp.value)}")
print_bp_entry(BadPixels.NOISE_OUT_OF_THRESHOLD)
print_bp_entry(BadPixels.OFFSET_NOISE_EVAL_ERROR)
print_bp_entry(BadPixels.WRONG_GAIN_VALUE)
def eval_bpidx(d):
mdn = np.nanmedian(d, axis=(0))[None, :, :]
std = np.nanstd(d, axis=(0))[None, :, :]
idx = (d > badpixel_threshold_sigma * std + mdn) | (
d < (-badpixel_threshold_sigma) * std + mdn
)
return idx
```
%% Cell type:code id:40c34cc5-fe93-4b83-bf39-f465f37c40b4 tags:
``` python
step_timer.start()
g_name = ["G0", "G1", "G2"]
for mod in karabo_da:
pdu = da_to_pdu[mod]
display(Markdown(f"### Badpixels for module {mod}:"))
badpixels_map[mod][
~np.isfinite(offset_map[mod])
] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
# Iteratively identify noisy pixels (OFFSET_NOISE_EVAL_ERROR & NOISE_OUT_OF_THRESHOLD)
# Due to the low number of strips (1280), outliers can significantly impact
# the badpixel calculation. We use multiple iterations to gradually refine
# the sigma threshold, as a single fixed threshold may not be suitable for
# different Gotthard2 detectors.
last_nbpix = 0
nmap = noise_map[mod].copy()
for itr in range(bpix_noise_nitrs):
badpixels_map[mod][
~np.isfinite(nmap)
] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
badpixels_map[mod][
eval_bpidx(nmap)
] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
nmap[badpixels_map[mod] > 0] = np.nan
recent_nbpix = np.count_nonzero(badpixels_map[mod])
if last_nbpix == recent_nbpix:
# stop iterating if no more badpixels difference
break
if not local_output:
for cell in [0, 1]:
fig, ax = plt.subplots(figsize=(10, 5))
for g_idx in [0, 1, 2]:
ax.plot(badpixels_map[mod][:, cell, g_idx], label=f"G{g_idx} Bad pixel map")
ax.set_xticks(np.arange(0, 1281, 80))
ax.set_xlabel("Stripes #")
ax.set_xlabel("BadPixels")
ax.set_title(f"BadPixels map - Cell {cell} - Module {mod} ({pdu})")
ax.set_ylim([0, 5])
ax.legend()
plt.show()
step_timer.done_step(f"Creating bad pixels constant.")
```
%% Cell type:code id:c8777cfe tags:
``` python
if not local_output:
for cons, cname in zip([offset_map, noise_map], ["Offset", "Noise"]):
for mod in karabo_da:
pdu = da_to_pdu[mod]
display(Markdown(f"### {cname} for module {mod}:"))
for cell in [0, 1]:
fig, ax = plt.subplots(figsize=(10, 5))
for g_idx in [0, 1, 2]:
ax.plot(cons[mod][:, cell, g_idx], label=f"G{g_idx} {cname} map")
ax.set_xlabel("Stripes #")
ax.set_xlabel(cname)
ax.set_title(f"{cname} map - Cell {cell} - Module {mod} ({pdu})")
ax.legend()
plt.show()
```
%% Cell type:code id:1c4eddf7-7d6e-49f4-8cbb-12d2bc496a8f tags:
``` python
step_timer.start()
for mod in karabo_da:
db_mod = da_to_pdu[mod]
constants = {
"Offset": offset_map[mod],
"Noise": noise_map[mod],
"BadPixelsDark": badpixels_map[mod],
}
md = None
for key, const_data in constants.items():
const = getattr(Constants.Gotthard2, key)()
const.data = const_data
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"• Exposure time: {exposure_time}\n"
f"• Exposure period: {exposure_period}\n"
f"• Acquisition rate: {acquisition_rate}\n"
f"• Single photon: {single_photon}\n"
f"• Creation time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\n"
)
step_timer.done_step("Injecting constants.")
```
%% Cell type:code id:98ca9486 tags:
``` python
# TODO: store old constants for comparison.
for mod in karabo_da:
mod_file = mod.replace("/", "-")
with open(f"{metadata_folder or out_folder}/module_metadata_{mod_file}.yml", "w") as fd:
yaml.safe_dump(
{
"module": mod,
"pdu": da_to_pdu[mod],
},
fd,
)
```
......
%% Cell type:markdown id:bed7bd15-21d9-4735-82c1-c27c1a5e3346 tags:
# Gotthard2 Offline Correction
Author: European XFEL Detector Group, Version: 1.0
Offline Correction for Gotthard2 Detector.
This notebook is able to correct 25um and 50um GH2 detectors using the same correction steps:
- Convert 12bit raw data into 10bit, offset subtraction, then multiply with gain constant.
| Correction | constants | boolean to enable/disable |
|------------|-------------|-----------------------------|
| 12bit to 10bit | `LUTGotthard2` | |
| Offset | `OffsetGotthard2`|`offset_correction`|
| Relative gain | `RelativeGainGotthard2` + `BadPixelsFFGotthard2` |`gain_correction`|
Beside the corrected data, a mask is stored using the badpixels constant of the same parameter conditions and time.
- `BadPixelsDarkGotthard2`
- `BadPixelsFFGotthard2`, if relative gain correction is requested.
The correction is done per sequence file. If all selected sequence files have no images to correct the notebook will fail.
The same result would be reached in case the needed dark calibration constants were not retrieved for all modules and `offset_correction` is True.
In case one of the gain constants were not retrieved `gain_correction` is switched to False and gain correction is disabled.
The `data` datasets stored in the RECEIVER source along with the corrected image (`adc`) and `mask` are:
- `gain`
- `bunchId`
- `memoryCell`
- `frameNumber`
- `timestamp`
- `trainId`
%% Cell type:code id:570322ed-f611-4fd1-b2ec-c12c13d55843 tags:
``` python
in_folder = "/gpfs/exfel/exp/DETLAB/202330/p900326/raw" # the folder to read data from, required
out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/gotthard2" # the folder to output to, required
metadata_folder = "" # Directory containing calibration_metadata.yml when run by xfel-calibrate
run = 20 # 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.
output_source_template = "{karabo_id}/CORR/RECEIVER:daqOutput" # Correction data source template. filled with karabo_id and correction receiver
ctrl_source_template = "{karabo_id_control}/DET/{control_template}" # template for control source name (filled with karabo_id_control)
karabo_id = "DETLAB_25UM_GH2" # karabo prefix of Gotthard-II devices
karabo_da = [""] # data aggregators
input_source_affixes = [""] # The affix to format into the input source template to be able to load the correct module's data.
input_source_template = "{karabo_id}/DET/RECEIVER{input_source_affixes}:daqOutput" # data source template used to read INSTRUMENT keys.
control_template = "CONTROL" # control template used to read CONTROL keys.
karabo_id_control = "" # Control karabo ID. Set to empty string to use the karabo-id
# Parameters for calibration database.
cal_db_interface = "tcp://max-exfl-cal001:8016#8025" # the database interface to use.
cal_db_timeout = 180000 # timeout on caldb requests.
creation_time = "" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC e.g. "2022-06-28 13:00:00"
# Parameters affecting corrected data.
offset_correction = True # apply offset correction. This can be disabled to only apply LUT or apply LUT and gain correction for non-linear differential results.
gain_correction = True # apply gain correction.
chunks_data = 1 # HDF chunk size for pixel data in number of frames.
# 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.
acquisition_rate = -1. # Detector acquisition rate (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.
reverse_second_module = -1 # Reverse 25um GH2 second module before interleaving. set to -1 to use value in raw file to reverse based on `CTRL/reverseSlaveReadOutMode`'s value.
sensor_type = "" # Specifies the GH2 sensor type: '50um', '25um', or "" to automatically decide from RUN data.
# 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 logging import warning
import h5py
import pasha as psh
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
import cal_tools.restful_config as rest_cfg
from cal_tools.calcat_interface import CalCatError, GOTTHARD2_CalibrationData
from cal_tools.files import DataFile
from cal_tools.gotthard2 import gotthard2algs, gotthard2lib
from cal_tools.step_timing import StepTimer
from cal_tools.tools import (
calcat_creation_time,
exit_notebook,
map_seq_files,
write_constants_fragment,
)
from XFELDetAna.plotting.heatmap import heatmapPlot
warnings.filterwarnings('ignore')
%matplotlib inline
```
%% Cell type:code id:fc893eff tags:
``` python
# Input parameters validation
gotthard2lib.validate_sensor_type(sensor_type)
```
%% 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
ctrl_src = ctrl_source_template.format(
karabo_id_control=karabo_id_control,
control_template=control_template,
)
# Run's creation time:
creation_time = calcat_creation_time(in_folder, run, creation_time)
print(f"Creation time: {creation_time}")
```
%% Cell type:code id:f9a8d1eb-ce6a-4ed0-abf4-4a6029734672 tags:
``` python
step_timer = StepTimer()
```
%% Cell type:code id:892172d8 tags:
``` python
run_dc = RunDirectory(run_folder)
# Read slow data
g2ctrl = gotthard2lib.Gotthard2Ctrl(run_dc=run_dc, ctrl_src=ctrl_src)
if bias_voltage == -1:
bias_voltage = g2ctrl.get_bias_voltage()
if exposure_time == -1:
exposure_time = g2ctrl.get_exposure_time()
if exposure_period == -1:
exposure_period = g2ctrl.get_exposure_period()
if acquisition_rate == -1:
acquisition_rate = g2ctrl.get_acquisition_rate()
if single_photon == -1:
single_photon = g2ctrl.get_single_photon()
if not sensor_type:
sensor_type = g2ctrl.get_sensor_type()
# Check if the GH2 detector type is "50um" and validate that only one data aggregator is provided
if sensor_type == "50um" and len(karabo_da) > 1:
raise ValueError(
f"karabo_da has more than one aggregator {karabo_da}, "
"but the GH2 detector type is specified as '50um', which expects only one aggregator."
)
gh2_detector = g2ctrl.get_det_type()
if reverse_second_module == -1 and gh2_detector == "25um":
if reverse_second_module == -1 and sensor_type == "25um":
reverse_second_module = not g2ctrl.second_module_reversed()
if reverse_second_module:
warning(
"Second module is not reversed. "
"Reversing second module before correction.")
print("Bias Voltage:", bias_voltage)
print("Exposure Time:", exposure_time)
print("Exposure Period:", exposure_period)
print("Acquisition Rate:", acquisition_rate)
print("Single Photon:", single_photon)
print(f"Processing {gh2_detector} Gotthard2.")
print(f"Processing {sensor_type} Gotthard2.")
```
%% Cell type:code id:21a8953a-8c76-475e-8f4f-b201cc25c159 tags:
``` python
# GH2 calibration data object.
g2_cal = GOTTHARD2_CalibrationData(
detector_name=karabo_id,
sensor_bias_voltage=bias_voltage,
exposure_time=exposure_time,
exposure_period=exposure_period,
acquisition_rate=acquisition_rate,
single_photon=single_photon,
event_at=creation_time,
client=rest_cfg.calibration_client(),
)
da_to_pdu = None
# Keep as long as it is essential to correct
# RAW data (FXE p003225) before the data mapping was added to CALCAT.
try: # in case local constants are used with old RAW data. This can be removed in the future.
da_to_pdu = g2_cal.mod_to_pdu
except CalCatError as e:
print(e)
db_modules = [None] * len(karabo_da)
if da_to_pdu:
if karabo_da == [""]:
karabo_da = sorted(da_to_pdu.keys())
else:
# Exclude non selected DA from processing.
karabo_da = [da for da in karabo_da if da in da_to_pdu]
db_modules = [da_to_pdu[da] for da in karabo_da]
print(f"Process modules: {db_modules} for run {run}")
da_to_source = gotthard2lib.map_da_to_source(
run_dc,
karabo_da,
input_source_template,
karabo_id,
input_source_affixes,
)
data_sources = list(da_to_source.values())
```
%% Cell type:code id:2551b923 tags:
``` python
# Check the available trains to correct.
total_trains = len(
RunDirectory(run_folder).select(data_sources, require_all=True).train_ids)
if total_trains:
print(f"Correcting {total_trains}.")
else:
raise ValueError(f"No trains to correct for run {run}.")
```
%% Cell type:markdown id:8c852392-bb19-4c40-b2ce-3b787538a92d tags:
### Retrieving calibration constants
%% Cell type:code id:5717d722 tags:
``` python
# Used for old FXE (p003225) runs before adding Gotthard2 to CALCAT
const_data = dict()
constant_names = ["LUTGotthard2", "OffsetGotthard2", "BadPixelsDarkGotthard2"]
if gain_correction:
constant_names += ["RelativeGainGotthard2", "BadPixelsFFGotthard2"]
g2_metadata = g2_cal.metadata(calibrations=constant_names)
# Display retrieved calibration constants timestamps
g2_cal.display_markdown_retrieved_constants(metadata=g2_metadata)
# Validate the constants availability and raise/warn correspondingly.
for mod, calibrations in g2_metadata.items():
dark_constants = {"LUTGotthard2"}
if offset_correction:
dark_constants |= {"OffsetGotthard2", "BadPixelsDarkGotthard2"}
missing_dark_constants = dark_constants - set(calibrations)
if missing_dark_constants:
karabo_da.remove(mod)
warning(f"Dark constants {missing_dark_constants} are not available to correct {mod}.") # noqa
missing_gain_constants = {
"BadPixelsFFGotthard2", "RelativeGainGotthard2"} - set(calibrations)
if gain_correction and missing_gain_constants:
warning(f"Gain constants {missing_gain_constants} are not retrieved for mod {mod}.")
if not karabo_da:
raise ValueError("Dark constants are not available for all modules.")
```
%% Cell type:code id:ac1cdec5 tags:
``` python
# Record constant details in YAML metadata.
write_constants_fragment(
out_folder=(metadata_folder or out_folder),
det_metadata=g2_metadata,
caldb_root=g2_cal.caldb_root)
# Load constants data for all constants.
const_data = g2_cal.ndarray_map(metadata=g2_metadata)
# Prepare constant arrays.
for mod in karabo_da:
# Create the mask array.
bpix = const_data[mod].get("BadPixelsDarkGotthard2")
if bpix is None:
bpix = np.zeros((1280, 2, 3), dtype=np.uint32)
if const_data[mod].get("BadPixelsFFGotthard2") is not None:
bpix |= const_data[mod]["BadPixelsFFGotthard2"]
const_data[mod]["Mask"] = bpix
# Prepare empty arrays for missing constants.
if const_data[mod].get("OffsetGotthard2") is None:
const_data[mod]["OffsetGotthard2"] = np.zeros(
(1280, 2, 3), dtype=np.float32)
if const_data[mod].get("RelativeGainGotthard2") is None:
const_data[mod]["RelativeGainGotthard2"] = np.ones(
(1280, 2, 3), dtype=np.float32)
const_data[mod]["RelativeGainGotthard2"] = const_data[mod]["RelativeGainGotthard2"].astype( # noqa
np.float32, copy=False) # Old gain constants are not float32.
```
%% Cell type:code id:2c7dd0bb tags:
``` python
file_da = list({kda.split('/')[0] for kda in karabo_da})
mapped_files, total_files = map_seq_files(
run_folder,
file_da,
sequences,
)
# This notebook doesn't account for processing more
# than one file data aggregator.
seq_files = mapped_files[file_da[0]]
if not len(seq_files):
raise IndexError(
"No sequence files available to correct for the selected sequences and karabo_da.")
print(f"Processing a total of {total_files} sequence files")
```
%% Cell type:code id:23fcf7f4-351a-4df7-8829-d8497d94fecc tags:
``` python
context = psh.ProcessContext(num_workers=23)
```
%% 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, const_data[mod]["LUTGotthard2"], data_corr[index, ...])
gotthard2algs.correct_train(
data_corr[index, ...],
mask[index, ...],
g,
const_data[mod]["OffsetGotthard2"].astype(np.float32), # PSI map is in f8
const_data[mod]["RelativeGainGotthard2"],
const_data[mod]["Mask"],
apply_offset=offset_correction,
apply_gain=gain_correction,
)
```
%% Cell type:code id:f88c1aa6-a735-4b72-adce-b30162f5daea tags:
``` python
corr_data_source = output_source_template.format(karabo_id=karabo_id)
empty_seq = 0
for raw_file in seq_files:
out_file = out_folder / raw_file.name.replace("RAW", "CORR")
# Select module INSTRUMENT sources and deselect empty trains.
dc = H5File(raw_file).select(data_sources, require_all=True)
n_trains = len(dc.train_ids)
if n_trains == 0:
empty_seq += 1
warning(f"Skipping correction. No trains to correct for this sequence file: {raw_file}.")
continue
else:
print(f"Correcting {n_trains} for {raw_file}.")
# Initialize GH2 data and gain arrays to store in corrected files.
if gh2_detector == "25um":
if sensor_type == "25um":
dshape_stored = (dc[data_sources[0], "data.adc"].shape[:2] + (1280 * 2,))
data_stored = np.zeros(dshape_stored, dtype=np.float32)
gain_stored = np.zeros(dshape_stored, dtype=np.uint8)
mask_stored = np.zeros(dshape_stored, dtype=np.uint32)
else:
data_stored, gain_stored, mask_stored = None, None, None
for i, mod in enumerate(karabo_da):
src = da_to_source[mod]
step_timer.start()
print(f"Correcting {src} for {raw_file}")
data = dc[src, "data.adc"].ndarray()
gain = dc[src, "data.gain"].ndarray()
# This is a safeguard, addressing the rare scenario of
# configuring the second module to not reverse the readout.
if reverse_second_module and i == 1:
data = np.flip(data, axis=-1)
gain = np.flip(gain, axis=-1)
step_timer.done_step("Preparing raw data")
dshape = data.shape
step_timer.start()
# 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(f"Correcting one receiver in one sequence file")
step_timer.start()
# Provided PSI gain map has 0 values. Set inf values to nan.
# TODO: This can maybe be removed after creating XFEL gain maps.?
data_corr[np.isinf(data_corr)] = np.nan
# Create CORR files and add corrected data sections.
image_counts = dc[src, "data.adc"].data_counts(labelled=False)
if gh2_detector == "25um":
if sensor_type == "25um":
data_stored[..., i::2] = data_corr
gain_stored[..., i::2] = gain
mask_stored[..., i::2] = mask
else: # "50um"
data_stored, gain_stored, mask_stored = data_corr, gain, mask
seq_file = dc.files[0] # FileAccess
with DataFile(out_file, "w") as ofile:
# Create INDEX datasets.
ofile.create_index(dc.train_ids, from_file=seq_file)
ofile.create_metadata(
like=dc,
sequence=seq_file.sequence,
instrument_channels=(f"{corr_data_source}/data",)
)
# Create Instrument section to later add corrected datasets.
outp_source = ofile.create_instrument_source(corr_data_source)
# Create count/first datasets at INDEX source.
outp_source.create_index(data=image_counts)
# Store uncorrected trainId in the corrected file.
outp_source.create_key(
f"data.trainId", data=dc.train_ids,
chunks=min(50, len(dc.train_ids))
)
# Create datasets with the available corrected data
for field_name, field_data in {
"adc": data_stored,
"gain": gain_stored,
}.items():
outp_source.create_key(
f"data.{field_name}", data=field_data,
chunks=((chunks_data,) + data_corr.shape[1:])
)
# For GH2 25um, the data of the second receiver is
# stored in the corrected file.
for field in ["bunchId", "memoryCell", "frameNumber", "timestamp"]:
outp_source.create_key(
f"data.{field}", data=dc[src, f"data.{field}"].ndarray(),
chunks=(chunks_data, data_corr.shape[1])
)
outp_source.create_compressed_key(f"data.mask", data=mask_stored)
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:
exit_notebook('Skipping plots as configured.')
if empty_seq == len(seq_files):
exit_notebook(f"All sequence files contain no data for correction.")
```
%% Cell type:code id:ff203f77-3811-46f3-bf7d-226d2dcab13f tags:
``` python
mod_dcs = {}
first_seq_raw = seq_files[0]
first_seq_corr = out_folder / first_seq_raw.name.replace("RAW", "CORR")
mod_dcs[corr_data_source] = {}
with H5File(first_seq_corr) as out_dc:
tid, mod_dcs[corr_data_source]["train_corr_data"] = next(
out_dc[corr_data_source, "data.adc"].trains()
)
mask = out_dc[corr_data_source, "data.mask"].train_from_id(tid)[1]
if gh2_detector == "25um":
if sensor_type == "25um":
mod_dcs[corr_data_source]["train_raw_data"] = np.zeros((data_corr.shape[1], 1280 * 2), dtype=np.float32)
mod_dcs[corr_data_source]["train_raw_gain"] = np.zeros((data_corr.shape[1], 1280 * 2), dtype=np.uint8)
for i, src in enumerate(data_sources):
with H5File(first_seq_raw) as in_dc:
train_dict = in_dc.train_from_id(tid)[1][src]
if gh2_detector == "25um":
if sensor_type == "25um":
if reverse_second_module and i == 1:
data = np.flip(train_dict["data.adc"], axis=-1)
gain = np.flip(train_dict["data.gain"], axis=-1)
else:
data = train_dict["data.adc"]
gain = train_dict["data.gain"]
mod_dcs[corr_data_source]["train_raw_data"][..., i::2] = data
mod_dcs[corr_data_source]["train_raw_gain"][..., i::2] = gain
else:
mod_dcs[corr_data_source]["train_raw_data"] = train_dict["data.adc"]
mod_dcs[corr_data_source]["train_raw_gain"] = train_dict["data.gain"]
```
%% Cell type:code id:bcb2c472 tags:
``` python
def plot_mean_data(data, title_format, y_label, fig_size=(10, 6)):
_, ax = plt.subplots(figsize=fig_size)
mean_data = np.mean(data, axis=0)
ax.plot(mean_data)
ax.set_title(title_format, fontsize=12)
ax.set_xlabel("Strip #", size=12)
ax.set_ylabel(y_label, size=12)
plt.xticks(fontsize=10)
plt.yticks(fontsize=10)
plt.show()
# Display markdown headers
display(Markdown("### Mean RAW and CORRECTED across pulses for one train:"))
display(Markdown(f"Train: {tid}"))
# Set title format based on gh2_detector
if gh2_detector == "50um":
# Set title format based on sensor_type
if sensor_type == "50um":
title_format = f"{{}} data for {karabo_da} ({db_modules})"
else:
title_format = f"Interleaved {{}} data for {karabo_da} ({db_modules})"
step_timer.start()
# Plot RAW data
raw_data = mod_dcs[corr_data_source]["train_raw_data"]
plot_mean_data(raw_data, title_format.format("RAW"), "12-bit ADC output")
# Plot CORRECTED data
corr_data = mod_dcs[corr_data_source]["train_corr_data"]
plot_mean_data(corr_data, title_format.format("CORRECTED"), "10-bit KeV. output")
# Plot CORRECTED MASKED data
corr_data_masked = corr_data.copy()
corr_data_masked[mask!=0] = np.nan
plot_mean_data(corr_data_masked, title_format.format("CORRECTED MASKED"), "10-bit KeV. output")
step_timer.done_step("Plotting mean data")
```
%% Cell type:code id:58a6a276 tags:
``` python
display(Markdown(f"### RAW and CORRECTED strips across pulses for train {tid}"))
step_timer.start()
for plt_data, dname in zip(
[raw_data, corr_data, corr_data_masked], ["RAW", "CORRECTED", "CORRECTED MASKED"]
):
fig, ax = plt.subplots(figsize=(15, 10))
plt.rcParams.update({"font.size": 14})
vmin, vmax = np.nanpercentile(plt_data, [1, 99])
cb_aspect = max(2, (vmax-vmin) / 25)
ax.tick_params(axis='both', which='major', labelsize=14)
heatmapPlot(
plt_data,
y_label="Pulses",
x_label="Strips",
title=title_format.format(dname),
use_axis=ax,
cb_aspect=cb_aspect,
vmin=vmin, vmax=vmax,
cmap='viridis',
)
pass
step_timer.done_step("Plotting RAW and CORRECTED data for one train")
```
%% 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:e2ecafe1 tags:
``` python
def plot_pulses(data, odd_pulse, even_pulse, title, y_label, fig_size=(10, 6)):
fig, ax = plt.subplots(figsize=fig_size)
ax.plot(data[odd_pulse], label=f"Odd Pulse {odd_pulse}")
if even_pulse:
ax.plot(data[even_pulse], label=f"Even Pulse {even_pulse}")
ax.set_title(title, fontsize=12)
ax.set_xlabel("Strip #", size=12)
ax.set_ylabel(y_label, size=12)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
ax.legend()
plt.show()
# Display markdown headers
display(Markdown("### RAW and CORRECTED even/odd pulses for one train:"))
display(Markdown(f"Train: {tid}"))
# Plot RAW data
plot_pulses(
raw_data, odd_pulse, even_pulse, title_format.format("RAW"), "12-bit ADC RAW")
# Plot CORRECTED data
plot_pulses(
corr_data, odd_pulse, even_pulse,
title_format.format("CORRECTED"), "10-bit KeV CORRECTED")
# Plot CORRECTED MASKED data
plot_pulses(
corr_data_masked, odd_pulse, even_pulse,
title_format.format("CORRECTED MASKED"), "10-bit KeV CORRECTED MASKED")
step_timer.done_step("Plotting RAW and CORRECTED odd/even pulses.")
```
......
%% Cell type:markdown id: tags:
# Gotthard2 Dark Summary
Author: European XFEL Detector Department, Version: 1.0
Summary for process dark constants and a comparison with previously injected constants with the same conditions.
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/DETLAB/202330/p900326/raw" # the folder to read data from, required
out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/gotthard2" # path to output to, required
metadata_folder = "" # Directory containing calibration_metadata.yml when run by xfel-calibrate.
run_high = 20 # run number for G0 dark run, required
run_med = 21 # run number for G1 dark run, required
run_low = 22 # run number for G2 dark run, required
# Parameters used to access raw data.
karabo_id = "DETLAB_25UM_GH2" # detector identifier.
karabo_da = ["DA01/1", "DA01/2"] # list of data aggregators, which corresponds to different JF modules. This is only needed for the detectors of one module.
control_template = "CONTROL" # control template used to read CONTROL keys.
ctrl_source_template = '{karabo_id}/DET/{control_template}' # 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
sensor_type = "" # Specifies the GH2 sensor type: '50um', '25um', or "" to automatically decide from RUN data.
# Parameters to be used for injecting dark calibration constants.
local_output = True # Boolean indicating that local constants were stored in the out_folder
# Skip the whole notebook if local_output is false in the preceding notebooks.
if not local_output:
from cal_tools.tools import exit_notebook
exit_notebook('No local constants saved. Skipping summary plots')
```
%% Cell type:code id: tags:
``` python
import warnings
from pathlib import Path
warnings.filterwarnings('ignore')
import h5py
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import yaml
from extra_data import RunDirectory
from IPython.display import Markdown, display
from cal_tools.gotthard2 import gotthard2lib
matplotlib.use("agg")
%matplotlib inline
from cal_tools.tools import CalibrationMetadata
```
%% Cell type:code id: tags:
``` python
# Input parameters validation
gotthard2lib.validate_sensor_type(sensor_type)
```
%% Cell type:code id: tags:
``` python
out_folder = Path(out_folder)
metadata = CalibrationMetadata(metadata_folder or out_folder)
```
%% Cell type:code id: tags:
``` python
if not karabo_id_control:
karabo_id_control = karabo_id
g2ctrl = gotthard2lib.Gotthard2Ctrl(
run_dc=RunDirectory(Path(in_folder) / f"r{run_high:04d}"),
ctrl_src=ctrl_source_template.format(
karabo_id_control=karabo_id_control,
control_template=control_template
))
gh2_detector = g2ctrl.get_det_type()
print(f"Processing {gh2_detector} Gotthard2.")
if not sensor_type:
sensor_type = g2ctrl.get_sensor_type()
print(f"Processing {sensor_type} Gotthard2.")
```
%% Cell type:code id: tags:
``` python
mod_mapping = dict()
for fn in Path(metadata_folder or out_folder).glob("module_metadata_*.yml"):
with fn.open("r") as fd:
fdict = yaml.safe_load(fd)
mod_mapping[fdict["module"].replace("-", "/")] = fdict["pdu"]
mod_mapping = dict(sorted(mod_mapping.items(), key=lambda item: item[0]))
```
%% Cell type:code id: tags:
``` python
plt_map = dict()
dark_constants = ["Offset", "Noise", "BadPixelsDark"]
if gh2_detector == "25um":
if sensor_type == "25um":
for cname in dark_constants:
plt_map[cname] = np.zeros(
(1280 * 2, 2, 3),
dtype=np.uint32 if cname == "BadPixelsDark" else np.float32
)
```
%% Cell type:code id: tags:
``` python
for cname in dark_constants:
for i, (mod, pdu) in enumerate(mod_mapping.items()):
module = mod.replace("-", "/")
with h5py.File(out_folder / f"const_{cname}_{pdu}.h5", 'r') as f:
if gh2_detector == "25um":
if sensor_type == "25um":
plt_map[cname][i::2] = f["data"][()]
else:
plt_map[cname] = f["data"][()]
```
%% Cell type:code id: tags:
``` python
if gh2_detector == "50um":
if sensor_type == "50um":
title = (
f"{{}} data for "
f"{[f'{mod}({pdu})' for mod, pdu in mod_mapping.items()]}"
)
else:
title = (
f"Interleaved {{}} data for "
f"{[f'{mod}({pdu})' for mod, pdu in mod_mapping.items()]}"
)
for cname in dark_constants:
display(Markdown(f"### {cname}"))
for cell in [0, 1]:
fig, ax = plt.subplots(figsize=(10, 5))
for g_idx in [0, 1, 2]:
ax.plot(plt_map[cname][:, cell, g_idx], label=f"G{g_idx} {cname} map")
ax.set_xticks(
np.arange(0, plt_map[cname].shape[0]+1, plt_map[cname].shape[0]//16)
)
ax.set_xlabel("Stripes #", fontsize=15)
ax.set_ylabel("ADU", fontsize=15)
ax.set_title(title.format(f"{cname} map - Cell {cell}"), fontsize=15)
if cname == "BadPixelsDark":
ax.set_ylim([0, 5])
ax.legend()
plt.show()
```
......
......@@ -61,7 +61,7 @@ install_requires = [
"dill==0.3.8",
"docutils==0.20.1",
"dynaconf==3.2.4",
"dynflatfield==1.0.0",
"dynflatfield==1.1.0",
"env_cache==0.1",
"extra_data==1.16.0",
"extra_geom==1.11.0",
......
......@@ -3,6 +3,16 @@ from pathlib import Path
from typing import List, Union
import extra_data
import numpy as np
def validate_sensor_type(sensor_type):
choices = ["25um", "50um"]
if sensor_type and sensor_type not in choices:
raise ValueError(
f"Unexpected Gotthard2 `sensor_type` given: '{sensor_type}'. "
f"Valid options are: {', '.join(choices)}."
)
def map_da_to_source(dc, das, source_template, karabo_id, affixes):
......@@ -63,16 +73,29 @@ class Gotthard2Ctrl():
return bool(
self.run_dc[self.ctrl_src, "singlePhoton"].as_single_value())
def get_det_type(self):
def get_sensor_type(self):
"""GH2 rxHostname is a vector of bytes objects.
GH2 25um has two host names unlike 50um which has one.
Returns:
str: return if its a 25um or 50um detector.
"""
hostname = self.run_dc.get_run_value(self.ctrl_src, "rxHostname")
return "25um" if hostname[1].decode("utf-8") else "50um"
hostnames = self.run_dc.get_run_value(self.ctrl_src, "rxHostname")
num_hosts = np.count_nonzero(hostnames != b'')
if num_hosts == 0:
raise ValueError(
"No hosts are present in RUN/rxHostname. "
"Unable to determine the sensor type for this data. ")
if num_hosts == 1:
return "50um"
elif num_hosts == 2:
return "25um"
else:
raise ValueError(
f"RUN/rxHostname has more host names ({num_hosts}) than expected. "
"Expected either 1 hostname for 50um or 2 host names for 25um."
)
def second_module_reversed(self):
"""Check if reverseSlaveReadOutMode is True."""
......