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

Merge branch 'feat/new_agipd_ff_localoutput' into 'master'

[AGIPD][FLAT FIELDS] Initial Implementation

See merge request detectors/pycalibration!380
parents 72eb7a35 19a1c95d
No related branches found
No related tags found
1 merge request!380[AGIPD][FLAT FIELDS] Initial Implementation
...@@ -624,7 +624,7 @@ class AgipdCorrections: ...@@ -624,7 +624,7 @@ class AgipdCorrections:
# slopeFF = slopeFFpix/avarege(slopeFFpix) # slopeFF = slopeFFpix/avarege(slopeFFpix)
# To apply them we have to / not * # To apply them we have to / not *
if self.corr_bools.get("xray_corr"): if self.corr_bools.get("xray_corr"):
data /= self.xray_cor[module_idx] data /= self.xray_cor[module_idx][cellid, ...]
# use sharedmem raw_data and t0_rgain # use sharedmem raw_data and t0_rgain
# after calculating it while offset correcting. # after calculating it while offset correcting.
...@@ -818,7 +818,7 @@ class AgipdCorrections: ...@@ -818,7 +818,7 @@ class AgipdCorrections:
uq, fidxv, cntsv = np.unique(trains, return_index=True, uq, fidxv, cntsv = np.unique(trains, return_index=True,
return_counts=True) return_counts=True)
# Validate calculated CORR INDEX contents by checking # Validate calculated CORR INDEX contents by checking
# difference between trainId stored in RAW data and trains from # difference between trainId stored in RAW data and trains from
train_diff = np.isin(np.array(infile["/INDEX/trainId"]), uq, train_diff = np.isin(np.array(infile["/INDEX/trainId"]), uq,
invert=True) invert=True)
...@@ -905,12 +905,12 @@ class AgipdCorrections: ...@@ -905,12 +905,12 @@ class AgipdCorrections:
exists of the current AGIPD instances. exists of the current AGIPD instances.
Relative gain is derived both from pulse capacitor as well as low Relative gain is derived both from pulse capacitor as well as low
intensity flat field data, information from flat field data is intensity flat field data, information from flat field data is
needed to 'calibrate' pulse capacitor data, if there is no needed to 'calibrate' pulse capacitor data, if there is no
available FF data, relative gain for High Gain stage is set to 1: available FF data, relative gain for High Gain stage is set to 1:
* Relative gain for High gain stage - from the FF data we get * Relative gain for High gain stage - from the FF data we get
the relative slopes of a given pixel and memory cells with the relative slopes of a given pixel and memory cells with
respect to all memory cells and all pixels in the module, respect to all memory cells and all pixels in the module,
Please note: Current slopesFF avaialble in calibibration Please note: Current slopesFF avaialble in calibibration
constants are created per pixel only, not per memory cell: constants are created per pixel only, not per memory cell:
...@@ -923,9 +923,9 @@ class AgipdCorrections: ...@@ -923,9 +923,9 @@ class AgipdCorrections:
between high and medium gain using slope information from between high and medium gain using slope information from
fits to the linear part of high and medium gain: fits to the linear part of high and medium gain:
rfpc_high_medium = m_h/m_m rfpc_high_medium = m_h/m_m
where m_h and m_m is the medium gain slope of given memory cells where m_h and m_m is the medium gain slope of given memory cells
and pixel and m_h is the high gain slope as above and pixel and m_h is the high gain slope as above
rel_gain_medium = rel_high_gain * rfpc_high_medium rel_gain_medium = rel_high_gain * rfpc_high_medium
...@@ -954,32 +954,16 @@ class AgipdCorrections: ...@@ -954,32 +954,16 @@ class AgipdCorrections:
if self.corr_bools.get("xray_corr"): if self.corr_bools.get("xray_corr"):
bpixels |= cons_data["BadPixelsFF"].astype(np.uint32)[..., :bpixels.shape[2], None] # noqa bpixels |= cons_data["BadPixelsFF"].astype(np.uint32)[..., :bpixels.shape[2], None] # noqa
slopesFF = np.squeeze(cons_data["SlopesFF"]) slopesFF = cons_data["SlopesFF"]
if len(slopesFF.shape) == 4: if len(slopesFF.shape) == 4:
slopesFF = slopesFF[..., 0] slopesFF = slopesFF[..., 0]
# Memory cell resolved xray_cor correction
# first 32 cells are known to behave differently so if we can avoid xray_cor = slopesFF # (128, 512, mem_cells)
# them
# when calculating the mean X-ray derived gain slope for each pixel
if slopesFF.shape[2] > 32:
xray_cor = np.nanmedian(
slopesFF[..., 32:min(96, self.max_cells)], axis=2)
elif slopesFF.shape[2] > 2:
xray_cor = np.nanmedian(
slopesFF[..., :min(96, self.max_cells)], axis=2)
else:
xray_cor = np.squeeze(slopesFF[..., 0])
# relative X-ray correction is normalized by the median # relative X-ray correction is normalized by the median
# of all pixels # of all pixels
# TODO: A check is required to know why it is again divided by xray_cor /= np.nanmedian(xray_cor)
# median. If we have relative slopes in the constants
# and (we have!)
# xray cor = (slopeFF/avarege_slopeFF)/avarege_slopeFF.
# It didn't not make sense and was removed.
# xray_cor /= np.nanmedian(xray_cor)
self.xray_cor[module_idx][...] = xray_cor.transpose()[...] self.xray_cor[module_idx][...] = xray_cor.transpose()[...]
...@@ -1041,13 +1025,19 @@ class AgipdCorrections: ...@@ -1041,13 +1025,19 @@ class AgipdCorrections:
# ration between HG and MG per pixel per mem cell used # ration between HG and MG per pixel per mem cell used
# for rel gain calculation # for rel gain calculation
frac_high_med_pix = pc_high_m / pc_med_m frac_high_med_pix = pc_high_m / pc_med_m
# avarage ration between HG and MG as a function of # avarage ration between HG and MG as a function of
# mem cell (needed for bls_stripes) # mem cell (needed for bls_stripes)
# TODO: Per pixel would be more optimal correction # TODO: Per pixel would be more optimal correction
frac_high_med = pc_high_med / pc_med_med frac_high_med = pc_high_med / pc_med_med
# calculate additional medium-gain offset # calculate additional medium-gain offset
md_additional_offset = pc_high_l - pc_med_l * pc_high_m / pc_med_m md_additional_offset = pc_high_l - pc_med_l * pc_high_m / pc_med_m
# Calculate relative gain. If FF constants are available,
# use them for high gain
# if not rel_gain is calculated using PC data only
# if self.corr_bools.get("xray_corr"):
# rel_gain[..., :self.max_cells, 0] /= xray_corr
# PC data should be 'calibrated with X-ray data, # PC data should be 'calibrated with X-ray data,
# if it is not done, it is better to use 1 instead of bias # if it is not done, it is better to use 1 instead of bias
# the results with PC arteffacts. # the results with PC arteffacts.
...@@ -1078,17 +1068,35 @@ class AgipdCorrections: ...@@ -1078,17 +1068,35 @@ class AgipdCorrections:
dname = device.device_name dname = device.device_name
cons_data = dict() cons_data = dict()
when = dict() when = dict()
for cname, mdata in const_yaml[dname].items(): for cname, mdata in const_yaml[dname].items():
when[cname] = mdata["creation-time"] when[cname] = mdata["creation-time"]
if when[cname]: if when[cname]:
cf = h5py.File(mdata["file-path"], "r") # This path is only used when testing new flat fields from
cons_data[cname] = np.copy(cf[f"{dname}/{cname}/0/data"]) # file during development: it takes ages to test using all
cf.close() # cells. Consequently, the shape needs to be fixed when less
# cells are used.
with h5py.File(mdata["file-path"], "r") as cf:
cons_data[cname] = np.copy(cf[f"{dname}/{cname}/0/data"])
shape = cons_data[cname].shape # (128, 512, mem_cells)
extra_dims = shape[:2] + (self.max_cells-shape[2], )
if extra_dims[-1] != 0 and cname == "BadPixelsFF":
extra_temp = np.zeros(extra_dims, dtype=np.int32)
cons_data[cname] = np.concatenate(
(cons_data[cname], extra_temp), axis=2)
print('An extra dimension was added to the constants '
'for the benefit of BadPixelsFF')
if extra_dims[-1] != 0 and cname == "SlopesFF":
extra_temp = np.ones(extra_dims, dtype=np.float32)
cons_data[cname] = np.concatenate(
(cons_data[cname], extra_temp), axis=2)
print('An extra dimension was added to the constants '
'for the benefit of SlopesFF')
else: else:
# Create empty constant using the list elements # Create empty constant using the list elements
cons_data[cname] = \ cons_data[cname] = getattr(np, mdata["file-path"][0])(mdata["file-path"][1]) # noqa
getattr(np, mdata["file-path"][0])(mdata["file-path"][1])
self.init_constants(cons_data, module_idx) self.init_constants(cons_data, module_idx)
return when return when
...@@ -1191,7 +1199,7 @@ class AgipdCorrections: ...@@ -1191,7 +1199,7 @@ class AgipdCorrections:
dtype='f4') dtype='f4')
self.mask[module_idx] = sharedmem.empty(constant_shape, dtype='i4') self.mask[module_idx] = sharedmem.empty(constant_shape, dtype='i4')
self.xray_cor[module_idx] = sharedmem.empty(constant_shape[2:], self.xray_cor[module_idx] = sharedmem.empty(constant_shape[1:],
dtype='f4') dtype='f4')
def allocate_images(self, shape, n_cores_files): def allocate_images(self, shape, n_cores_files):
......
...@@ -43,7 +43,7 @@ def assemble_constant_dict(corr_bools, pc_bools, memory_cells, bias_voltage, ...@@ -43,7 +43,7 @@ def assemble_constant_dict(corr_bools, pc_bools, memory_cells, bias_voltage,
const_dict = { const_dict = {
"Offset": ["zeros", (128, 512, memory_cells, 3), darkcond], "Offset": ["zeros", (128, 512, memory_cells, 3), darkcond],
"Noise": ["zeros", (128, 512, memory_cells, 3), darkcond], "Noise": ["zeros", (128, 512, memory_cells, 3), darkcond],
"ThresholdsDark": ["ones", (128, 512, memory_cells, 2), darkcond], "ThresholdsDark": ["ones", (128, 512, memory_cells, 5), darkcond],
"BadPixelsDark": ["zeros", (128, 512, memory_cells, 3), darkcond], "BadPixelsDark": ["zeros", (128, 512, memory_cells, 3), darkcond],
} }
......
from typing import Any, Dict, List, Optional, Tuple
from iminuit import Minuit
import numpy as np
from cal_tools.enums import BadPixelsFF
def any_in(mask: np.ndarray, bits: int) -> bool:
return mask.astype(np.uint) & bits > 0
def gaussian(x: np.ndarray, norm: int = 1, mean: int = 0, sigma: int = 1) -> float: # noqa
"""
Return value of Gaussian function
:param x: Argument (float of 1D array) of Gaussian function
:param norm: Normalization of Gaussian function
:param mean: Mean parameter
:param sigma: Sigma parameter
:return: Value of gaussian function.
"""
return norm * np.exp(-1 / 2 * ((x - mean) / sigma) ** 2) / (sigma * np.sqrt(2 * np.pi)) # noqa
def gaussian_sum(x: np.ndarray, ng: int = 4, *p: Tuple[Any]) -> float:
"""Sum of ng Gaussian functions
:param x: Argument (float of 1D array) of the function
:param ng: Number of Gaussian functions
:param p: List of parameters (norm1,mean1,sigma1,norm2,mean2,sigma2,...)
"""
r = 0.
for i in range(ng):
r += gaussian(x, *p[i * 3:(i + 1) * 3])
return r
def get_statistical_parameters(x: np.ndarray,
y: np.ndarray,
x_range: List) -> Tuple[np.uint64, np.float64, np.float64, np.ndarray]: # noqa
"""Return statistical parameters of selected part of a histogram.
:param x: Center of bins of the histogram
:param y: Value of bins of the histogram
:param x_range: x range to be considered
:return: Sum of histogram, Mean value, Standard deviation,
List of selected bins
"""
# TODO: Check if wq.median works better than mean
sel = (x >= x_range[0]) & (x < x_range[1])
h_sum = np.sum(y[sel])
h_norm = y[sel] / h_sum
h_mean = np.sum(h_norm * x[sel])
h_sqr = (x[sel] - h_mean) ** 2
h_std = np.sqrt(np.sum(h_norm * h_sqr))
return h_sum, h_mean, h_std, sel
def get_starting_parameters(xe: np.ndarray,
ye: np.ndarray,
limits: np.ndarray,
n_peaks: int = 3,
f_lim: int = 2) -> Tuple[Dict[str, Any], List[Tuple]]: # noqa
"""
Estimate starting parameters for Gaussian fit of several peaks.
:param xe: Center of bins of the histogram
:param ye: Value of bins of the histogram
:param limits: Position of each peak ((left1, right1),
(left2, right2), ...) to be considered.
:param n_peaks: Number of peaks
:param f_lim: Limits in units of standard deviation to consider
"""
parameters = {}
shapes = []
for peak in range(n_peaks):
n, m, rms, idx = get_statistical_parameters(xe, ye, limits[peak])
limits2 = [m - f_lim * rms, m + f_lim * rms]
n, m, rms, idx = get_statistical_parameters(xe, ye, limits2)
shapes.append((n, m, rms, idx))
parameters.update({f'g{peak}sigma': rms,
f'g{peak}n': float(n),
f'g{peak}mean': m})
return parameters, shapes
def fit_n_peaks(x: np.ndarray,
y: np.ndarray,
pars: Dict[str, Any],
x_range: Tuple[float, float],
do_minos: Optional[bool] = False,
n_peaks: Optional[int] = 4,
fix_d01: Optional[bool] = True) -> Minuit:
"""
Fit histogram with n-Gaussian function.
:param x: Center of bins of the histogram
:param y: Value of bins of the histogram
:param pars: Dictionary of initial parameters for fitting
:param x_range: x Range to be considered for the fitting
:param do_minos: Run Minos if True
:param n_peaks: Number of Gaussian peaks to fit (min 2, max 4)
:param fix_d01: Fix position of peaks to the distance between noise and
first photon peak.
:return: minuit object
"""
sel = (x >= x_range[0]) & (x < x_range[1])
# Square of bin errors
yrr2 = np.copy(y[sel])
yrr2[yrr2 == 0] = 1 # bins with zero events have error=1
if fix_d01:
pars['fix_g2mean'] = True
pars['fix_g3mean'] = True
if n_peaks < 4:
pars['g3n'] = 0
pars['fix_g3n'] = True
pars['g3sigma'] = 1
pars['fix_g3sigma'] = True
pars['fix_g3mean'] = True
if n_peaks < 3:
pars['g2n'] = 0
pars['fix_g2n'] = True
pars['g2sigma'] = 1
pars['fix_g2sigma'] = True
pars['fix_g2mean'] = True
def chi2_f(g0n, g0mean, g0sigma,
g1n, g1mean, g1sigma,
g2n, g2mean, g2sigma,
g3n, g3mean, g3sigma, ):
d01 = (g1mean - g0mean)
if 'fix_g2mean' in pars and pars['fix_g2mean']:
g2mean = g0mean + d01 * 2
if 'fix_g3mean' in pars and pars['fix_g3mean']:
g3mean = g0mean + d01 * 3
if g3n == 0:
n_peaks = 3
elif g2n == 0:
n_peaks = 2
else:
n_peaks = 4
yt = gaussian_sum(x[sel], n_peaks,
g0n, g0mean, g0sigma,
g1n, g1mean, g1sigma,
g2n, g2mean, g2sigma,
g3n, g3mean, g3sigma)
return np.nansum((yt - y[sel]) ** 2 / yrr2)
minuit = Minuit(chi2_f, **pars, pedantic=False)
minuit.migrad()
if do_minos:
if minuit.get_fmin().is_valid:
minuit.minos()
return minuit
def set_par_limits(pars: Dict[str, Any],
peak_range: np.ndarray,
peak_norm_range: np.ndarray,
peak_width_range: np.ndarray,
n_peaks: Optional[int] = 4):
"""
Set limits on initial fit parameters based on given values
:param pars: Dictionary of initial fit parameters
:param peak_range: Limits on peak positions
:param peak_norm_range: Limits on normalization of Gaussian peaks
:param peak_width_range: Limits on width of Gaussian peaks
:param n_peaks: Number of Gaussian peaks
"""
for peak in range(n_peaks):
pars.update({f'limit_g{peak}n': peak_norm_range[peak],
f'limit_g{peak}mean': peak_range[peak],
f'limit_g{peak}sigma': peak_width_range[peak],
})
def get_mask(fit_summary: Dict[str, Any],
peak_lim: List,
d0_lim: List,
chi2_lim: List,
peak_width_lim: np.array) -> int:
"""
Calculate Bad pixels mask based on fit results and given limits
:param fit_summary: Dictionary of the fit output from Multi-Gaussian fit
:param peak_lim: Limits on noise peak position
:param d0_lim: Limits on distance between noise and first photon peak
:param chi2_lim: Limits on reduced chi^2 value
:param peak_width_lim: Limits on noise peak width
:return: Bad pixel mask
"""
if not fit_summary['is_valid']:
return BadPixelsFF.FIT_FAILED.value
m0 = fit_summary['g0mean']
s0 = fit_summary['g0sigma']
s1 = fit_summary['g1sigma']
s2 = fit_summary['g2sigma']
chi2_ndof = fit_summary['chi2_ndof']
d01 = fit_summary['g1mean'] - m0
mask = 0
if not fit_summary['is_valid']:
mask |= BadPixelsFF.FIT_FAILED.value
if not fit_summary['has_accurate_covar']:
mask |= BadPixelsFF.ACCURATE_COVAR.value
if not peak_lim[0] <= m0 <= peak_lim[1]:
mask |= BadPixelsFF.NOISE_PEAK_THRESHOLD.value
if not d0_lim[0] <= d01 <= d0_lim[1]:
mask |= BadPixelsFF.GAIN_THRESHOLD.value
if not chi2_lim[0] <= chi2_ndof <= chi2_lim[1]:
mask |= BadPixelsFF.CHI2_THRESHOLD.value
width_lim = peak_width_lim[0] * s0
inside_s1 = width_lim[0] <= s1 <= width_lim[1]
width_lim = peak_width_lim[1] * s0
inside_s2 = width_lim[0] <= s2 <= width_lim[1]
if not inside_s1 and inside_s2:
mask |= BadPixelsFF.PEAK_WIDTH_THRESHOLD.value
return mask
...@@ -28,6 +28,21 @@ class BadPixels(Enum): ...@@ -28,6 +28,21 @@ class BadPixels(Enum):
NON_LIN_RESPONSE_REGION = 0b100000000000000000000 # bit 21 NON_LIN_RESPONSE_REGION = 0b100000000000000000000 # bit 21
class BadPixelsFF(Enum):
""" The SLopesFF Bad Pixel Encoding
"""
FIT_FAILED = 0b000000000000000000001 # bit 1
CHI2_THRESHOLD = 0b000000000000000000010 # bit 2
NOISE_PEAK_THRESHOLD = 0b000000000000000000100 # bit 3
GAIN_THRESHOLD = 0b000000000000000001000 # bit 4
PEAK_WIDTH_THRESHOLD = 0b000000000000000010000 # bit 5
ACCURATE_COVAR = 0b000000000000000100000 # bit 6
BAD_DARK = 0b000000000000001000000 # bit 6
NO_ENTRY = 0b000000000000010000000 # bit 7
GAIN_DEVIATION = 0b000000000000100000000 # bit 8
class SnowResolution(Enum): class SnowResolution(Enum):
""" An Enum specifying how to resolve snowy pixels """ An Enum specifying how to resolve snowy pixels
""" """
......
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# AGIPD Offline Correction # # AGIPD Offline Correction #
Author: European XFEL Detector Group, Version: 2.0 Author: European XFEL Detector Group, Version: 2.0
Offline Calibration for the AGIPD Detector Offline Calibration for the AGIPD Detector
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
in_folder = "/gpfs/exfel/exp/HED/202031/p900174/raw" # the folder to read data from, required in_folder = "/gpfs/exfel/exp/HED/202031/p900174/raw" # the folder to read data from, required
out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/hibef_agipd2" # the folder to output to, required out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/hibef_agipd2" # the folder to output to, required
sequences = [-1] # sequences to correct, set to -1 for all, range allowed sequences = [-1] # sequences to correct, set to -1 for all, range allowed
modules = [-1] # modules to correct, set to -1 for all, range allowed modules = [-1] # modules to correct, set to -1 for all, range allowed
run = 155 # runs to process, required run = 155 # runs to process, required
karabo_id = "HED_DET_AGIPD500K2G" # karabo karabo_id karabo_id = "HED_DET_AGIPD500K2G" # karabo karabo_id
karabo_da = ['-1'] # a list of data aggregators names, Default [-1] for selecting all data aggregators karabo_da = ['-1'] # a list of data aggregators names, Default [-1] for selecting all data aggregators
receiver_id = "{}CH0" # inset for receiver devices receiver_id = "{}CH0" # inset for receiver devices
path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data
h5path = 'INSTRUMENT/{}/DET/{}:xtdf/' # path in the HDF5 file to images h5path = 'INSTRUMENT/{}/DET/{}:xtdf/' # path in the HDF5 file to images
h5path_idx = 'INDEX/{}/DET/{}:xtdf/' # path in the HDF5 file to images h5path_idx = 'INDEX/{}/DET/{}:xtdf/' # path in the HDF5 file to images
h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP' # path to control information h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP' # path to control information
karabo_id_control = "HED_EXP_AGIPD500K2G" # karabo-id for control device karabo_id_control = "HED_EXP_AGIPD500K2G" # karabo-id for control device
karabo_da_control = 'AGIPD500K2G00' # karabo DA for control infromation karabo_da_control = 'AGIPD500K2G00' # karabo DA for control infromation
slopes_ff_from_files = "" # Path to locally stored SlopesFF and BadPixelsFF constants
use_dir_creation_date = True # use the creation data of the input dir for database queries use_dir_creation_date = True # use the creation data of the input dir for database queries
cal_db_interface = "tcp://max-exfl016:8015#8045" # the database interface to use cal_db_interface = "tcp://max-exfl016:8015#8045" # the database interface to use
cal_db_timeout = 30000 # in milli seconds cal_db_timeout = 30000 # in milli seconds
creation_date_offset = "00:00:00" # add an offset to creation date, e.g. to get different constants creation_date_offset = "00:00:00" # add an offset to creation date, e.g. to get different constants
max_cells = 0 # number of memory cells used, set to 0 to automatically infer max_cells = 0 # number of memory cells used, set to 0 to automatically infer
bias_voltage = 300 # Bias voltage bias_voltage = 300 # Bias voltage
acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine
gain_setting = 0.1 # the gain setting, use 0.1 to try to auto-determine gain_setting = 0.1 # the gain setting, use 0.1 to try to auto-determine
photon_energy = 9.2 # photon energy in keV photon_energy = 9.2 # photon energy in keV
overwrite = True # set to True if existing data should be overwritten overwrite = True # set to True if existing data should be overwritten
max_pulses = [0, 500, 1] # range list [st, end, step] of maximum pulse indices within a train. 3 allowed maximum list input elements. max_pulses = [0, 500, 1] # range list [st, end, step] of maximum pulse indices within a train. 3 allowed maximum list input elements.
mem_cells_db = 0 # set to a value different than 0 to use this value for DB queries mem_cells_db = 0 # set to a value different than 0 to use this value for DB queries
cell_id_preview = 1 # cell Id used for preview in single-shot plots cell_id_preview = 1 # cell Id used for preview in single-shot plots
# Correction parameters # Correction parameters
blc_noise_threshold = 5000 # above this mean signal intensity now baseline correction via noise is attempted blc_noise_threshold = 5000 # above this mean signal intensity now baseline correction via noise is attempted
cm_dark_fraction = 0.66 # threshold for fraction of empty pixels to consider module enough dark to perform CM correction cm_dark_fraction = 0.66 # threshold for fraction of empty pixels to consider module enough dark to perform CM correction
cm_dark_range = [-50.,30] # range for signal value ADU for pixel to be consider as a dark pixel cm_dark_range = [-50.,30] # range for signal value ADU for pixel to be consider as a dark pixel
cm_n_itr = 4 # number of iterations for common mode correction cm_n_itr = 4 # number of iterations for common mode correction
hg_hard_threshold = 1000 # threshold to force medium gain offset subtracted pixel to high gain hg_hard_threshold = 1000 # threshold to force medium gain offset subtracted pixel to high gain
mg_hard_threshold = 1000 # threshold to force medium gain offset subtracted pixel from low to medium gain mg_hard_threshold = 1000 # threshold to force medium gain offset subtracted pixel from low to medium gain
noisy_adc_threshold = 0.25 # threshold to mask complete adc noisy_adc_threshold = 0.25 # threshold to mask complete adc
# Correction Booleans # Correction Booleans
only_offset = False # Apply only Offset correction. if False, Offset is applied by Default. if True, Offset is only applied. only_offset = False # Apply only Offset correction. if False, Offset is applied by Default. if True, Offset is only applied.
rel_gain = False # do relative gain correction based on PC data rel_gain = False # do relative gain correction based on PC data
xray_gain = False # do relative gain correction based on xray data xray_gain = False # do relative gain correction based on xray data
blc_noise = False # if set, baseline correction via noise peak location is attempted blc_noise = False # if set, baseline correction via noise peak location is attempted
blc_stripes = False # if set, baseline corrected via stripes blc_stripes = False # if set, baseline corrected via stripes
blc_hmatch = False # if set, base line correction via histogram matching is attempted blc_hmatch = False # if set, base line correction via histogram matching is attempted
match_asics = False # if set, inner ASIC borders are matched to the same signal level match_asics = False # if set, inner ASIC borders are matched to the same signal level
adjust_mg_baseline = False # adjust medium gain baseline to match highest high gain value adjust_mg_baseline = False # adjust medium gain baseline to match highest high gain value
zero_nans = False # set NaN values in corrected data to 0 zero_nans = False # set NaN values in corrected data to 0
zero_orange = False # set to 0 very negative and very large values in corrected data zero_orange = False # set to 0 very negative and very large values in corrected data
blc_set_min = False # Shift to 0 negative medium gain pixels after offset corr blc_set_min = False # Shift to 0 negative medium gain pixels after offset corr
corr_asic_diag = False # if set, diagonal drop offs on ASICs are correted corr_asic_diag = False # if set, diagonal drop offs on ASICs are correted
force_hg_if_below = False # set high gain if mg offset subtracted value is below hg_hard_threshold force_hg_if_below = False # set high gain if mg offset subtracted value is below hg_hard_threshold
force_mg_if_below = False # set medium gain if mg offset subtracted value is below mg_hard_threshold force_mg_if_below = False # set medium gain if mg offset subtracted value is below mg_hard_threshold
mask_noisy_adc = False # Mask entire ADC if they are noise above a relative threshold mask_noisy_adc = False # Mask entire ADC if they are noise above a relative threshold
common_mode = False # Common mode correction common_mode = False # Common mode correction
melt_snow = False # Identify (and optionally interpolate) 'snowy' pixels melt_snow = False # Identify (and optionally interpolate) 'snowy' pixels
mask_zero_std = False # Mask pixels with zero standard deviation across train mask_zero_std = False # Mask pixels with zero standard deviation across train
low_medium_gap = False # 5 sigma separation in thresholding between low and medium gain low_medium_gap = False # 5 sigma separation in thresholding between low and medium gain
# Paralellization parameters # Paralellization parameters
chunk_size = 1000 # Size of chunk for image-weise correction chunk_size = 1000 # Size of chunk for image-weise correction
chunk_size_idim = 1 # chunking size of imaging dimension, adjust if user software is sensitive to this. chunk_size_idim = 1 # chunking size of imaging dimension, adjust if user software is sensitive to this.
n_cores_correct = 16 # Number of chunks to be processed in parallel n_cores_correct = 16 # Number of chunks to be processed in parallel
n_cores_files = 4 # Number of files to be processed in parallel n_cores_files = 4 # Number of files to be processed in parallel
sequences_per_node = 2 # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel sequences_per_node = 2 # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel
def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da): def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da):
from xfel_calibrate.calibrate import balance_sequences as bs from xfel_calibrate.calibrate import balance_sequences as bs
return bs(in_folder, run, sequences, sequences_per_node, karabo_da) return bs(in_folder, run, sequences, sequences_per_node, karabo_da)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import copy import copy
from datetime import timedelta from datetime import timedelta
from dateutil import parser from dateutil import parser
import gc import gc
import glob import glob
import itertools import itertools
from IPython.display import HTML, display, Markdown, Latex from IPython.display import HTML, display, Markdown, Latex
import math import math
from multiprocessing import Pool from multiprocessing import Pool
import os import os
import re import re
import sys import sys
import traceback import traceback
from time import time, sleep, perf_counter from time import time, sleep, perf_counter
import tabulate import tabulate
import warnings import warnings
warnings.filterwarnings('ignore') warnings.filterwarnings('ignore')
import yaml import yaml
from extra_geom import AGIPD_1MGeometry, AGIPD_500K2GGeometry from extra_geom import AGIPD_1MGeometry, AGIPD_500K2GGeometry
from extra_data import RunDirectory, stack_detector_data from extra_data import RunDirectory, stack_detector_data
from iCalibrationDB import Detectors from iCalibrationDB import Detectors
from mpl_toolkits.mplot3d import Axes3D from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import LinearLocator, FormatStrFormatter from matplotlib.ticker import LinearLocator, FormatStrFormatter
from matplotlib.colors import LogNorm from matplotlib.colors import LogNorm
from matplotlib import cm as colormap from matplotlib import cm as colormap
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib import matplotlib
matplotlib.use("agg") matplotlib.use("agg")
%matplotlib inline %matplotlib inline
import numpy as np import numpy as np
import seaborn as sns import seaborn as sns
sns.set() sns.set()
sns.set_context("paper", font_scale=1.4) sns.set_context("paper", font_scale=1.4)
sns.set_style("ticks") sns.set_style("ticks")
from cal_tools.agipdlib import (AgipdCorrections, get_acq_rate, from cal_tools.agipdlib import (AgipdCorrections, get_acq_rate,
get_gain_setting, get_num_cells) get_gain_setting, get_num_cells)
from cal_tools.cython import agipdalgs as calgs from cal_tools.cython import agipdalgs as calgs
from cal_tools.ana_tools import get_range from cal_tools.ana_tools import get_range
from cal_tools.enums import BadPixels from cal_tools.enums import BadPixels
from cal_tools.tools import get_dir_creation_date, map_modules_from_folder from cal_tools.tools import get_dir_creation_date, map_modules_from_folder
from cal_tools.step_timing import StepTimer from cal_tools.step_timing import StepTimer
import seaborn as sns import seaborn as sns
sns.set() sns.set()
sns.set_context("paper", font_scale=1.4) sns.set_context("paper", font_scale=1.4)
sns.set_style("ticks") sns.set_style("ticks")
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Evaluated parameters ## ## Evaluated parameters ##
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Fill dictionaries comprising bools and arguments for correction and data analysis # Fill dictionaries comprising bools and arguments for correction and data analysis
# Here the herarichy and dependability for correction booleans are defined # Here the herarichy and dependability for correction booleans are defined
corr_bools = {} corr_bools = {}
# offset is at the bottom of AGIPD correction pyramid. # offset is at the bottom of AGIPD correction pyramid.
corr_bools["only_offset"] = only_offset corr_bools["only_offset"] = only_offset
# Dont apply any corrections if only_offset is requested # Dont apply any corrections if only_offset is requested
if not only_offset: if not only_offset:
corr_bools["adjust_mg_baseline"] = adjust_mg_baseline corr_bools["adjust_mg_baseline"] = adjust_mg_baseline
corr_bools["rel_gain"] = rel_gain corr_bools["rel_gain"] = rel_gain
corr_bools["xray_corr"] = xray_gain corr_bools["xray_corr"] = xray_gain
corr_bools["blc_noise"] = blc_noise corr_bools["blc_noise"] = blc_noise
corr_bools["blc_stripes"] = blc_stripes corr_bools["blc_stripes"] = blc_stripes
corr_bools["blc_hmatch"] = blc_hmatch corr_bools["blc_hmatch"] = blc_hmatch
corr_bools["blc_set_min"] = blc_set_min corr_bools["blc_set_min"] = blc_set_min
corr_bools["match_asics"] = match_asics corr_bools["match_asics"] = match_asics
corr_bools["corr_asic_diag"] = corr_asic_diag corr_bools["corr_asic_diag"] = corr_asic_diag
corr_bools["zero_nans"] = zero_nans corr_bools["zero_nans"] = zero_nans
corr_bools["zero_orange"] = zero_orange corr_bools["zero_orange"] = zero_orange
corr_bools["mask_noisy_adc"] = mask_noisy_adc corr_bools["mask_noisy_adc"] = mask_noisy_adc
corr_bools["force_hg_if_below"] = force_hg_if_below corr_bools["force_hg_if_below"] = force_hg_if_below
corr_bools["force_mg_if_below"] = force_mg_if_below corr_bools["force_mg_if_below"] = force_mg_if_below
corr_bools["common_mode"] = common_mode corr_bools["common_mode"] = common_mode
corr_bools["melt_snow"] = melt_snow corr_bools["melt_snow"] = melt_snow
corr_bools["mask_zero_std"] = mask_zero_std corr_bools["mask_zero_std"] = mask_zero_std
corr_bools["low_medium_gap"] = low_medium_gap corr_bools["low_medium_gap"] = low_medium_gap
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if in_folder[-1] == "/": if in_folder[-1] == "/":
in_folder = in_folder[:-1] in_folder = in_folder[:-1]
if sequences[0] == -1: if sequences[0] == -1:
sequences = None sequences = None
control_fname = f'{in_folder}/r{run:04d}/RAW-R{run:04d}-{karabo_da_control}-S00000.h5' control_fname = f'{in_folder}/r{run:04d}/RAW-R{run:04d}-{karabo_da_control}-S00000.h5'
h5path_ctrl = h5path_ctrl.format(karabo_id_control) h5path_ctrl = h5path_ctrl.format(karabo_id_control)
h5path = h5path.format(karabo_id, receiver_id) h5path = h5path.format(karabo_id, receiver_id)
h5path_idx = h5path_idx.format(karabo_id, receiver_id) h5path_idx = h5path_idx.format(karabo_id, receiver_id)
print(f'Path to control file {control_fname}') print(f'Path to control file {control_fname}')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Create output folder # Create output folder
os.makedirs(out_folder, exist_ok=overwrite) os.makedirs(out_folder, exist_ok=overwrite)
# Evaluate detector instance for mapping # Evaluate detector instance for mapping
instrument = karabo_id.split("_")[0] instrument = karabo_id.split("_")[0]
if instrument == "SPB": if instrument == "SPB":
dinstance = "AGIPD1M1" dinstance = "AGIPD1M1"
nmods = 16 nmods = 16
elif instrument == "MID": elif instrument == "MID":
dinstance = "AGIPD1M2" dinstance = "AGIPD1M2"
nmods = 16 nmods = 16
# TODO: Remove DETLAB # TODO: Remove DETLAB
elif instrument == "HED" or instrument == "DETLAB": elif instrument == "HED" or instrument == "DETLAB":
dinstance = "AGIPD500K" dinstance = "AGIPD500K"
nmods = 8 nmods = 8
# Evaluate requested modules # Evaluate requested modules
if karabo_da[0] == '-1': if karabo_da[0] == '-1':
if modules[0] == -1: if modules[0] == -1:
modules = list(range(nmods)) modules = list(range(nmods))
karabo_da = ["AGIPD{:02d}".format(i) for i in modules] karabo_da = ["AGIPD{:02d}".format(i) for i in modules]
else: else:
modules = [int(x[-2:]) for x in karabo_da] modules = [int(x[-2:]) for x in karabo_da]
def mod_name(modno): def mod_name(modno):
return f"Q{modno // 4 + 1}M{modno % 4 + 1}" return f"Q{modno // 4 + 1}M{modno % 4 + 1}"
print("Process modules: ", ', '.join( print("Process modules: ", ', '.join(
[mod_name(x) for x in modules])) [mod_name(x) for x in modules]))
print(f"Detector in use is {karabo_id}") print(f"Detector in use is {karabo_id}")
print(f"Instrument {instrument}") print(f"Instrument {instrument}")
print(f"Detector instance {dinstance}") print(f"Detector instance {dinstance}")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Display Information about the selected pulses indices for correction. # Display Information about the selected pulses indices for correction.
pulses_lst = list(range(*max_pulses)) if not (len(max_pulses)==1 and max_pulses[0]==0) else max_pulses pulses_lst = list(range(*max_pulses)) if not (len(max_pulses)==1 and max_pulses[0]==0) else max_pulses
try: try:
if len(pulses_lst) > 1: if len(pulses_lst) > 1:
print("A range of {} pulse indices is selected: from {} to {} with a step of {}" print("A range of {} pulse indices is selected: from {} to {} with a step of {}"
.format(len(pulses_lst), pulses_lst[0] , pulses_lst[-1] + (pulses_lst[1] - pulses_lst[0]), .format(len(pulses_lst), pulses_lst[0] , pulses_lst[-1] + (pulses_lst[1] - pulses_lst[0]),
pulses_lst[1] - pulses_lst[0])) pulses_lst[1] - pulses_lst[0]))
else: else:
print("one pulse is selected: a pulse of idx {}".format(pulses_lst[0])) print("one pulse is selected: a pulse of idx {}".format(pulses_lst[0]))
except Exception as e: except Exception as e:
raise ValueError('max_pulses input Error: {}'.format(e)) raise ValueError('max_pulses input Error: {}'.format(e))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# set everything up filewise # set everything up filewise
mmf = map_modules_from_folder(in_folder, run, path_template, mmf = map_modules_from_folder(in_folder, run, path_template,
karabo_da, sequences) karabo_da, sequences)
mapped_files, mod_ids, total_sequences, sequences_qm, _ = mmf mapped_files, mod_ids, total_sequences, sequences_qm, _ = mmf
file_list = [] file_list = []
# ToDo: Split table over pages # ToDo: Split table over pages
print(f"Processing a total of {total_sequences} sequence files in chunks of {n_cores_files}") print(f"Processing a total of {total_sequences} sequence files in chunks of {n_cores_files}")
table = [] table = []
ti = 0 ti = 0
for k, files in mapped_files.items(): for k, files in mapped_files.items():
i = 0 i = 0
for f in list(files.queue): for f in list(files.queue):
file_list.append(f) file_list.append(f)
if i == 0: if i == 0:
table.append((ti, k, i, f)) table.append((ti, k, i, f))
else: else:
table.append((ti, "", i, f)) table.append((ti, "", i, f))
i += 1 i += 1
ti += 1 ti += 1
md = display(Latex(tabulate.tabulate(table, tablefmt='latex', md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
headers=["#", "module", "# module", "file"]))) headers=["#", "module", "# module", "file"])))
file_list = sorted(file_list, key=lambda name: name[-10:]) file_list = sorted(file_list, key=lambda name: name[-10:])
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
filename = file_list[0] filename = file_list[0]
channel = int(re.findall(r".*-AGIPD([0-9]+)-.*", filename)[0]) channel = int(re.findall(r".*-AGIPD([0-9]+)-.*", filename)[0])
# Evaluate number of memory cells # Evaluate number of memory cells
mem_cells = get_num_cells(filename, karabo_id, channel) mem_cells = get_num_cells(filename, karabo_id, channel)
if mem_cells is None: if mem_cells is None:
raise ValueError(f"No raw images found in {filename}") raise ValueError(f"No raw images found in {filename}")
mem_cells_db = mem_cells if mem_cells_db == 0 else mem_cells_db mem_cells_db = mem_cells if mem_cells_db == 0 else mem_cells_db
max_cells = mem_cells if max_cells == 0 else max_cells max_cells = mem_cells if max_cells == 0 else max_cells
# Evaluate aquisition rate # Evaluate aquisition rate
if acq_rate == 0: if acq_rate == 0:
acq_rate = get_acq_rate((filename, karabo_id, channel)) acq_rate = get_acq_rate((filename, karabo_id, channel))
print(f"Maximum memory cells to calibrate: {max_cells}") print(f"Maximum memory cells to calibrate: {max_cells}")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Evaluate creation time # Evaluate creation time
creation_time = None creation_time = None
if use_dir_creation_date: if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run) creation_time = get_dir_creation_date(in_folder, run)
offset = parser.parse(creation_date_offset) offset = parser.parse(creation_date_offset)
delta = timedelta(hours=offset.hour, delta = timedelta(hours=offset.hour,
minutes=offset.minute, seconds=offset.second) minutes=offset.minute, seconds=offset.second)
creation_time += delta creation_time += delta
# Evaluate gain setting # Evaluate gain setting
if gain_setting == 0.1: if gain_setting == 0.1:
if creation_time.replace(tzinfo=None) < parser.parse('2020-01-31'): if creation_time.replace(tzinfo=None) < parser.parse('2020-01-31'):
print("Set gain-setting to None for runs taken before 2020-01-31") print("Set gain-setting to None for runs taken before 2020-01-31")
gain_setting = None gain_setting = None
else: else:
try: try:
gain_setting = get_gain_setting(control_fname, h5path_ctrl) gain_setting = get_gain_setting(control_fname, h5path_ctrl)
except Exception as e: except Exception as e:
print(f'ERROR: while reading gain setting from: \n{control_fname}') print(f'ERROR: while reading gain setting from: \n{control_fname}')
print(e) print(e)
print("Set gain setting to 0") print("Set gain setting to 0")
gain_setting = 0 gain_setting = 0
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
print(f"Using {creation_time} as creation time") print(f"Using {creation_time} as creation time")
print(f"Operating conditions are:\n• Bias voltage: {bias_voltage}\n• Memory cells: {mem_cells_db}\n" print(f"Operating conditions are:\n• Bias voltage: {bias_voltage}\n• Memory cells: {mem_cells_db}\n"
f"• Acquisition rate: {acq_rate}\n• Gain setting: {gain_setting}\n• Photon Energy: {photon_energy}\n") f"• Acquisition rate: {acq_rate}\n• Gain setting: {gain_setting}\n• Photon Energy: {photon_energy}\n")
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Data processing ## ## Data processing ##
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
agipd_corr = AgipdCorrections(max_cells, max_pulses, agipd_corr = AgipdCorrections(max_cells, max_pulses,
h5_data_path=h5path, h5_data_path=h5path,
h5_index_path=h5path_idx, h5_index_path=h5path_idx,
corr_bools=corr_bools) corr_bools=corr_bools)
agipd_corr.baseline_corr_noise_threshold = -blc_noise_threshold agipd_corr.baseline_corr_noise_threshold = -blc_noise_threshold
agipd_corr.hg_hard_threshold = hg_hard_threshold agipd_corr.hg_hard_threshold = hg_hard_threshold
agipd_corr.mg_hard_threshold = mg_hard_threshold agipd_corr.mg_hard_threshold = mg_hard_threshold
agipd_corr.cm_dark_min = cm_dark_range[0] agipd_corr.cm_dark_min = cm_dark_range[0]
agipd_corr.cm_dark_max = cm_dark_range[1] agipd_corr.cm_dark_max = cm_dark_range[1]
agipd_corr.cm_dark_fraction = cm_dark_fraction agipd_corr.cm_dark_fraction = cm_dark_fraction
agipd_corr.cm_n_itr = cm_n_itr agipd_corr.cm_n_itr = cm_n_itr
agipd_corr.noisy_adc_threshold = noisy_adc_threshold agipd_corr.noisy_adc_threshold = noisy_adc_threshold
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Retrieve calibration constants to RAM # Retrieve calibration constants to RAM
agipd_corr.allocate_constants(modules, (3, mem_cells_db, 512, 128)) agipd_corr.allocate_constants(modules, (3, mem_cells_db, 512, 128))
const_yaml = None const_yaml = None
if os.path.isfile(f'{out_folder}/retrieved_constants.yml'): if os.path.isfile(f'{out_folder}/retrieved_constants.yml'):
with open(f'{out_folder}/retrieved_constants.yml', "r") as f: with open(f'{out_folder}/retrieved_constants.yml', "r") as f:
const_yaml = yaml.safe_load(f.read()) const_yaml = yaml.safe_load(f.read())
# retrive constants # retrive constants
def retrieve_constants(mod): def retrieve_constants(mod):
""" """
Retrieve calibration constants and load them to shared memory Retrieve calibration constants and load them to shared memory
Metadata for constants is taken from yml file or retrieved from the DB Metadata for constants is taken from yml file or retrieved from the DB
""" """
device = getattr(getattr(Detectors, dinstance), mod_name(mod)) device = getattr(getattr(Detectors, dinstance), mod_name(mod))
err = '' err = ''
try: try:
# check if there is a yaml file in out_folder that has the device constants. # check if there is a yaml file in out_folder that has the device constants.
if const_yaml and device.device_name in const_yaml: if const_yaml and device.device_name in const_yaml:
when = agipd_corr.initialize_from_yaml(const_yaml, mod, device) when = agipd_corr.initialize_from_yaml(const_yaml, mod, device)
else: else:
when = agipd_corr.initialize_from_db(cal_db_interface, creation_time, mem_cells_db, bias_voltage, when = agipd_corr.initialize_from_db(cal_db_interface, creation_time, mem_cells_db, bias_voltage,
photon_energy, gain_setting, acq_rate, mod, device, False) photon_energy, gain_setting, acq_rate, mod, device, False)
except Exception as e: except Exception as e:
err = f"Error: {e}\nError traceback: {traceback.format_exc()}" err = f"Error: {e}\nError traceback: {traceback.format_exc()}"
when = None when = None
return err, mod, when, device.device_name return err, mod, when, device.device_name
ts = perf_counter() ts = perf_counter()
with Pool(processes=len(modules)) as pool: with Pool(processes=len(modules)) as pool:
const_out = pool.map(retrieve_constants, modules) const_out = pool.map(retrieve_constants, modules)
print(f"Constants were loaded in {perf_counter()-ts:.01f}s") print(f"Constants were loaded in {perf_counter()-ts:.01f}s")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# allocate memory for images and hists # allocate memory for images and hists
n_images_max = max_cells*256 n_images_max = max_cells*256
data_shape = (n_images_max, 512, 128) data_shape = (n_images_max, 512, 128)
agipd_corr.allocate_images(data_shape, n_cores_files) agipd_corr.allocate_images(data_shape, n_cores_files)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def batches(l, batch_size): def batches(l, batch_size):
"""Group a list into batches of (up to) batch_size elements""" """Group a list into batches of (up to) batch_size elements"""
start = 0 start = 0
while start < len(l): while start < len(l):
yield l[start:start + batch_size] yield l[start:start + batch_size]
start += batch_size start += batch_size
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def imagewise_chunks(img_counts): def imagewise_chunks(img_counts):
"""Break up the loaded data into chunks of up to chunk_size """Break up the loaded data into chunks of up to chunk_size
Yields (file data slot, start index, stop index) Yields (file data slot, start index, stop index)
""" """
for i_proc, n_img in enumerate(img_counts): for i_proc, n_img in enumerate(img_counts):
n_chunks = math.ceil(n_img / chunk_size) n_chunks = math.ceil(n_img / chunk_size)
for i in range(n_chunks): for i in range(n_chunks):
yield i_proc, i * n_img // n_chunks, (i+1) * n_img // n_chunks yield i_proc, i * n_img // n_chunks, (i+1) * n_img // n_chunks
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
step_timer = StepTimer() step_timer = StepTimer()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
with Pool() as pool: with Pool() as pool:
for file_batch in batches(file_list, n_cores_files): for file_batch in batches(file_list, n_cores_files):
# TODO: Move some printed output to logging or similar # TODO: Move some printed output to logging or similar
print(f"Processing next {len(file_batch)} files:") print(f"Processing next {len(file_batch)} files:")
for file_name in file_batch: for file_name in file_batch:
print(" ", file_name) print(" ", file_name)
step_timer.start() step_timer.start()
img_counts = pool.starmap(agipd_corr.read_file, enumerate(file_batch)) img_counts = pool.starmap(agipd_corr.read_file, enumerate(file_batch))
step_timer.done_step('Loading data from files') step_timer.done_step('Loading data from files')
# Evaluate zero-data-std mask # Evaluate zero-data-std mask
pool.starmap(agipd_corr.mask_zero_std, itertools.product( pool.starmap(agipd_corr.mask_zero_std, itertools.product(
range(len(file_batch)), np.array_split(np.arange(agipd_corr.max_cells), n_cores_correct) range(len(file_batch)), np.array_split(np.arange(agipd_corr.max_cells), n_cores_correct)
)) ))
step_timer.done_step('Mask 0 std') step_timer.done_step('Mask 0 std')
# Perform offset image-wise correction # Perform offset image-wise correction
pool.starmap(agipd_corr.offset_correction, imagewise_chunks(img_counts)) pool.starmap(agipd_corr.offset_correction, imagewise_chunks(img_counts))
step_timer.done_step("Offset correction") step_timer.done_step("Offset correction")
if blc_noise or blc_stripes or blc_hmatch: if blc_noise or blc_stripes or blc_hmatch:
# Perform image-wise correction # Perform image-wise correction
pool.starmap(agipd_corr.baseline_correction, imagewise_chunks(img_counts)) pool.starmap(agipd_corr.baseline_correction, imagewise_chunks(img_counts))
step_timer.done_step("Base-line shift correction") step_timer.done_step("Base-line shift correction")
if common_mode: if common_mode:
# Perform cross-file correction parallel over asics # Perform cross-file correction parallel over asics
pool.starmap(agipd_corr.cm_correction, itertools.product( pool.starmap(agipd_corr.cm_correction, itertools.product(
range(len(file_batch)), range(16) # 16 ASICs per module range(len(file_batch)), range(16) # 16 ASICs per module
)) ))
step_timer.done_step("Common-mode correction") step_timer.done_step("Common-mode correction")
# Perform image-wise correction # Perform image-wise correction
pool.starmap(agipd_corr.gain_correction, imagewise_chunks(img_counts)) pool.starmap(agipd_corr.gain_correction, imagewise_chunks(img_counts))
step_timer.done_step("Image-wise correction") step_timer.done_step("Image-wise correction")
# Save corrected data # Save corrected data
pool.starmap(agipd_corr.write_file, [ pool.starmap(agipd_corr.write_file, [
(i_proc, file_name, os.path.join(out_folder, os.path.basename(file_name).replace("RAW", "CORR"))) (i_proc, file_name, os.path.join(out_folder, os.path.basename(file_name).replace("RAW", "CORR")))
for i_proc, file_name in enumerate(file_batch) for i_proc, file_name in enumerate(file_batch)
]) ])
step_timer.done_step("Save") step_timer.done_step("Save")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
print(f"Correction of {len(file_list)} files is finished") print(f"Correction of {len(file_list)} files is finished")
print(f"Total processing time {step_timer.timespan():.01f} s") print(f"Total processing time {step_timer.timespan():.01f} s")
print(f"Timing summary per batch of {n_cores_files} files:") print(f"Timing summary per batch of {n_cores_files} files:")
step_timer.print_summary() step_timer.print_summary()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# if there is a yml file that means a leading notebook got processed # if there is a yml file that means a leading notebook got processed
# and the reporting would be generated from it. # and the reporting would be generated from it.
fst_print = True fst_print = True
to_store = [] to_store = []
line = [] line = []
for i, (error, modno, when, mod_dev) in enumerate(const_out): for i, (error, modno, when, mod_dev) in enumerate(const_out):
qm = mod_name(modno) qm = mod_name(modno)
# expose errors while applying correction # expose errors while applying correction
if error: if error:
print("Error: {}".format(error) ) print("Error: {}".format(error) )
if not const_yaml or mod_dev not in const_yaml: if not const_yaml or mod_dev not in const_yaml:
if fst_print: if fst_print:
print("Constants are retrieved with creation time: ") print("Constants are retrieved with creation time: ")
fst_print = False fst_print = False
line = [qm] line = [qm]
# If correction is crashed # If correction is crashed
if not error: if not error:
print(f"{qm}:") print(f"{qm}:")
for key, item in when.items(): for key, item in when.items():
if hasattr(item, 'strftime'): if hasattr(item, 'strftime'):
item = item.strftime('%y-%m-%d %H:%M') item = item.strftime('%y-%m-%d %H:%M')
when[key] = item when[key] = item
print('{:.<12s}'.format(key), item) print('{:.<12s}'.format(key), item)
# Store few time stamps if exists # Store few time stamps if exists
# Add NA to keep array structure # Add NA to keep array structure
for key in ['Offset', 'SlopesPC', 'SlopesFF']: for key in ['Offset', 'SlopesPC', 'SlopesFF']:
if when and key in when and when[key]: if when and key in when and when[key]:
line.append(when[key]) line.append(when[key])
else: else:
if error is not None: if error is not None:
line.append('Err') line.append('Err')
else: else:
line.append('NA') line.append('NA')
if len(line) > 0: if len(line) > 0:
to_store.append(line) to_store.append(line)
seq = sequences[0] if sequences else 0 seq = sequences[0] if sequences else 0
if len(to_store) > 0: if len(to_store) > 0:
with open(f"{out_folder}/retrieved_constants_s{seq}.yml","w") as fyml: with open(f"{out_folder}/retrieved_constants_s{seq}.yml","w") as fyml:
yaml.safe_dump({"time-summary": {f"S{seq}":to_store}}, fyml) yaml.safe_dump({"time-summary": {f"S{seq}":to_store}}, fyml)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def do_3d_plot(data, edges, x_axis, y_axis): def do_3d_plot(data, edges, x_axis, y_axis):
fig = plt.figure(figsize=(10, 10)) fig = plt.figure(figsize=(10, 10))
ax = fig.gca(projection='3d') ax = fig.gca(projection='3d')
# Make data. # Make data.
X = edges[0][:-1] X = edges[0][:-1]
Y = edges[1][:-1] Y = edges[1][:-1]
X, Y = np.meshgrid(X, Y) X, Y = np.meshgrid(X, Y)
Z = data.T Z = data.T
# Plot the surface. # Plot the surface.
surf = ax.plot_surface(X, Y, Z, cmap=colormap.coolwarm, surf = ax.plot_surface(X, Y, Z, cmap=colormap.coolwarm,
linewidth=0, antialiased=False) linewidth=0, antialiased=False)
ax.set_xlabel(x_axis) ax.set_xlabel(x_axis)
ax.set_ylabel(y_axis) ax.set_ylabel(y_axis)
ax.set_zlabel("Counts") ax.set_zlabel("Counts")
def do_2d_plot(data, edges, y_axis, x_axis): def do_2d_plot(data, edges, y_axis, x_axis):
fig = plt.figure(figsize=(10, 10)) fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
extent = [np.min(edges[1]), np.max(edges[1]), extent = [np.min(edges[1]), np.max(edges[1]),
np.min(edges[0]), np.max(edges[0])] np.min(edges[0]), np.max(edges[0])]
im = ax.imshow(data[::-1, :], extent=extent, aspect="auto", im = ax.imshow(data[::-1, :], extent=extent, aspect="auto",
norm=LogNorm(vmin=1, vmax=max(10, np.max(data)))) norm=LogNorm(vmin=1, vmax=max(10, np.max(data))))
ax.set_xlabel(x_axis) ax.set_xlabel(x_axis)
ax.set_ylabel(y_axis) ax.set_ylabel(y_axis)
cb = fig.colorbar(im) cb = fig.colorbar(im)
cb.set_label("Counts") cb.set_label("Counts")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def get_trains_data(run_folder, source, include, tid=None, path='*/DET/*', modules=16, fillvalue=np.nan): def get_trains_data(run_folder, source, include, tid=None, path='*/DET/*', modules=16, fillvalue=np.nan):
""" """
Load single train for all module Load single train for all module
:param run_folder: Path to folder with data :param run_folder: Path to folder with data
:param source: Data source to be loaded :param source: Data source to be loaded
:param include: Inset of file name to be considered :param include: Inset of file name to be considered
:param tid: Train Id to be loaded. First train is considered if None is given :param tid: Train Id to be loaded. First train is considered if None is given
:param path: Path to find image data inside h5 file :param path: Path to find image data inside h5 file
""" """
run_data = RunDirectory(run_folder, include) run_data = RunDirectory(run_folder, include)
if tid: if tid:
tid, data = run_data.select('*/DET/*', source).train_from_id(tid) tid, data = run_data.select('*/DET/*', source).train_from_id(tid)
return tid, stack_detector_data(train=data, data=source, fillvalue=fillvalue, modules=modules) return tid, stack_detector_data(train=data, data=source, fillvalue=fillvalue, modules=modules)
else: else:
for tid, data in run_data.select('*/DET/*', source).trains(require_all=True): for tid, data in run_data.select('*/DET/*', source).trains(require_all=True):
return tid, stack_detector_data(train=data, data=source, fillvalue=fillvalue, modules=modules) return tid, stack_detector_data(train=data, data=source, fillvalue=fillvalue, modules=modules)
return None, None return None, None
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if dinstance == "AGIPD500K": if dinstance == "AGIPD500K":
geom = AGIPD_500K2GGeometry.from_origin() geom = AGIPD_500K2GGeometry.from_origin()
else: else:
geom = AGIPD_1MGeometry.from_quad_positions(quad_pos=[ geom = AGIPD_1MGeometry.from_quad_positions(quad_pos=[
(-525, 625), (-525, 625),
(-550, -10), (-550, -10),
(520, -160), (520, -160),
(542.5, 475), (542.5, 475),
]) ])
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
include = '*S00000*' if sequences is None else f'*S{sequences[0]:05d}*' include = '*S00000*' if sequences is None else f'*S{sequences[0]:05d}*'
tid, corrected = get_trains_data(f'{out_folder}/', 'image.data', include, modules=nmods) tid, corrected = get_trains_data(f'{out_folder}/', 'image.data', include, modules=nmods)
_, gains = get_trains_data(f'{out_folder}/', 'image.gain', include, tid, modules=nmods) _, gains = get_trains_data(f'{out_folder}/', 'image.gain', include, tid, modules=nmods)
_, mask = get_trains_data(f'{out_folder}/', 'image.mask', include, tid, modules=nmods) _, mask = get_trains_data(f'{out_folder}/', 'image.mask', include, tid, modules=nmods)
_, blshift = get_trains_data(f'{out_folder}/', 'image.blShift', include, tid, modules=nmods) _, blshift = get_trains_data(f'{out_folder}/', 'image.blShift', include, tid, modules=nmods)
_, cellId = get_trains_data(f'{out_folder}/', 'image.cellId', include, tid, modules=nmods) _, cellId = get_trains_data(f'{out_folder}/', 'image.cellId', include, tid, modules=nmods)
_, pulseId = get_trains_data(f'{out_folder}/', 'image.pulseId', include, tid, _, pulseId = get_trains_data(f'{out_folder}/', 'image.pulseId', include, tid,
modules=nmods, fillvalue=0) modules=nmods, fillvalue=0)
_, raw = get_trains_data(f'{in_folder}/r{run:04d}/', 'image.data', include, tid, modules=nmods) _, raw = get_trains_data(f'{in_folder}/r{run:04d}/', 'image.data', include, tid, modules=nmods)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
display(Markdown(f'## Preview and statistics for {gains.shape[0]} images of the train {tid} ##\n')) display(Markdown(f'## Preview and statistics for {gains.shape[0]} images of the train {tid} ##\n'))
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Signal vs. Analogue Gain ### ### Signal vs. Analogue Gain ###
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
hist, bins_x, bins_y = calgs.histogram2d(raw[:,0,...].flatten().astype(np.float32), hist, bins_x, bins_y = calgs.histogram2d(raw[:,0,...].flatten().astype(np.float32),
raw[:,1,...].flatten().astype(np.float32), raw[:,1,...].flatten().astype(np.float32),
bins=(100, 100), bins=(100, 100),
range=[[4000, 8192], [4000, 8192]]) range=[[4000, 8192], [4000, 8192]])
do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Analogue gain (ADU)") do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Analogue gain (ADU)")
do_3d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Analogue gain (ADU)") do_3d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Analogue gain (ADU)")
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Signal vs. Digitized Gain ### ### Signal vs. Digitized Gain ###
The following plot shows plots signal vs. digitized gain The following plot shows plots signal vs. digitized gain
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
hist, bins_x, bins_y = calgs.histogram2d(corrected.flatten().astype(np.float32), hist, bins_x, bins_y = calgs.histogram2d(corrected.flatten().astype(np.float32),
gains.flatten().astype(np.float32), bins=(100, 3), gains.flatten().astype(np.float32), bins=(100, 3),
range=[[-50, 8192], [0, 3]]) range=[[-50, 8192], [0, 3]])
do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Gain bit value") do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Gain bit value")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
print(f"Gain statistics in %") print(f"Gain statistics in %")
table = [[f'{gains[gains==0].size/gains.size*100:.02f}', table = [[f'{gains[gains==0].size/gains.size*100:.02f}',
f'{gains[gains==1].size/gains.size*100:.03f}', f'{gains[gains==1].size/gains.size*100:.03f}',
f'{gains[gains==2].size/gains.size*100:.03f}']] f'{gains[gains==2].size/gains.size*100:.03f}']]
md = display(Latex(tabulate.tabulate(table, tablefmt='latex', md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
headers=["High", "Medium", "Low"]))) headers=["High", "Medium", "Low"])))
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Intensity per Pulse ### ### Intensity per Pulse ###
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
pulse_range = [np.min(pulseId[pulseId>=0]), np.max(pulseId[pulseId>=0])] pulse_range = [np.min(pulseId[pulseId>=0]), np.max(pulseId[pulseId>=0])]
mean_data = np.nanmean(corrected, axis=(2, 3)) mean_data = np.nanmean(corrected, axis=(2, 3))
hist, bins_x, bins_y = calgs.histogram2d(mean_data.flatten().astype(np.float32), hist, bins_x, bins_y = calgs.histogram2d(mean_data.flatten().astype(np.float32),
pulseId.flatten().astype(np.float32), pulseId.flatten().astype(np.float32),
bins=(100, int(pulse_range[1])), bins=(100, int(pulse_range[1])),
range=[[-50, 1000], pulse_range]) range=[[-50, 1000], pulse_range])
do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id") do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id")
do_3d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id") do_3d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id")
hist, bins_x, bins_y = calgs.histogram2d(mean_data.flatten().astype(np.float32), hist, bins_x, bins_y = calgs.histogram2d(mean_data.flatten().astype(np.float32),
pulseId.flatten().astype(np.float32), pulseId.flatten().astype(np.float32),
bins=(100, int(pulse_range[1])), bins=(100, int(pulse_range[1])),
range=[[-50, 200000], pulse_range]) range=[[-50, 200000], pulse_range])
do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id") do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id")
do_3d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id") do_3d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id")
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Baseline shift ### ### Baseline shift ###
Estimated base-line shift with respect to the total ADU counts of corrected image. Estimated base-line shift with respect to the total ADU counts of corrected image.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
h = ax.hist(blshift.flatten(), bins=100, log=True) h = ax.hist(blshift.flatten(), bins=100, log=True)
_ = plt.xlabel('Baseline shift [ADU]') _ = plt.xlabel('Baseline shift [ADU]')
_ = plt.ylabel('Counts') _ = plt.ylabel('Counts')
_ = ax.grid() _ = ax.grid()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(10, 10)) fig = plt.figure(figsize=(10, 10))
corrected_ave = np.nansum(corrected, axis=(2, 3)) corrected_ave = np.nansum(corrected, axis=(2, 3))
plt.scatter(corrected_ave.flatten()/10**6, blshift.flatten(), s=0.9) plt.scatter(corrected_ave.flatten()/10**6, blshift.flatten(), s=0.9)
plt.xlim(-1, 1000) plt.xlim(-1, 1000)
plt.grid() plt.grid()
plt.xlabel('Illuminated corrected [MADU] ') plt.xlabel('Illuminated corrected [MADU] ')
_ = plt.ylabel('Estimated baseline shift [ADU]') _ = plt.ylabel('Estimated baseline shift [ADU]')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
display(Markdown('### Raw preview ###\n')) display(Markdown('### Raw preview ###\n'))
display(Markdown(f'Mean over images of the RAW data\n')) display(Markdown(f'Mean over images of the RAW data\n'))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
data = np.mean(raw[:, 0, ...], axis=0) data = np.mean(raw[:, 0, ...], axis=0)
vmin, vmax = get_range(data, 5) vmin, vmax = get_range(data, 5)
ax = geom.plot_data_fast(data, ax=ax, cmap="jet", vmin=vmin, vmax=vmax) ax = geom.plot_data_fast(data, ax=ax, cmap="jet", vmin=vmin, vmax=vmax)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
display(Markdown(f'Single shot of the RAW data from cell {np.max(cellId[cell_id_preview])} \n')) display(Markdown(f'Single shot of the RAW data from cell {np.max(cellId[cell_id_preview])} \n'))
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
vmin, vmax = get_range(raw[cell_id_preview, 0, ...], 5) vmin, vmax = get_range(raw[cell_id_preview, 0, ...], 5)
ax = geom.plot_data_fast(raw[cell_id_preview, 0, ...], ax=ax, cmap="jet", vmin=vmin, vmax=vmax) ax = geom.plot_data_fast(raw[cell_id_preview, 0, ...], ax=ax, cmap="jet", vmin=vmin, vmax=vmax)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
display(Markdown('### Corrected preview ###\n')) display(Markdown('### Corrected preview ###\n'))
display(Markdown(f'A single shot image from cell {np.max(cellId[cell_id_preview])} \n')) display(Markdown(f'A single shot image from cell {np.max(cellId[cell_id_preview])} \n'))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
vmin, vmax = get_range(corrected[cell_id_preview], 7, -50) vmin, vmax = get_range(corrected[cell_id_preview], 7, -50)
vmin = - 50 vmin = - 50
ax = geom.plot_data_fast(corrected[cell_id_preview], ax=ax, cmap="jet", vmin=vmin, vmax=vmax) ax = geom.plot_data_fast(corrected[cell_id_preview], ax=ax, cmap="jet", vmin=vmin, vmax=vmax)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
vmin, vmax = get_range(corrected[cell_id_preview], 5, -50) vmin, vmax = get_range(corrected[cell_id_preview], 5, -50)
nbins = np.int((vmax + 50) / 2) nbins = np.int((vmax + 50) / 2)
h = ax.hist(corrected[cell_id_preview].flatten(), h = ax.hist(corrected[cell_id_preview].flatten(),
bins=nbins, range=(-50, vmax), bins=nbins, range=(-50, vmax),
histtype='stepfilled', log=True) histtype='stepfilled', log=True)
_ = plt.xlabel('[ADU]') _ = plt.xlabel('[ADU]')
_ = plt.ylabel('Counts') _ = plt.ylabel('Counts')
_ = ax.grid() _ = ax.grid()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
display(Markdown('### Mean CORRECTED Preview ###\n')) display(Markdown('### Mean CORRECTED Preview ###\n'))
display(Markdown(f'A mean across one train \n')) display(Markdown(f'A mean across one train \n'))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
data = np.mean(corrected, axis=0) data = np.mean(corrected, axis=0)
vmin, vmax = get_range(data, 7) vmin, vmax = get_range(data, 7)
ax = geom.plot_data_fast(data, ax=ax, cmap="jet", vmin=-50, vmax=vmax) ax = geom.plot_data_fast(data, ax=ax, cmap="jet", vmin=-50, vmax=vmax)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
vmin, vmax = get_range(corrected, 10, -100) vmin, vmax = get_range(corrected, 10, -100)
vmax = np.nanmax(corrected) vmax = np.nanmax(corrected)
if vmax > 50000: if vmax > 50000:
vmax=50000 vmax=50000
nbins = np.int((vmax + 100) / 5) nbins = np.int((vmax + 100) / 5)
h = ax.hist(corrected.flatten(), bins=nbins, h = ax.hist(corrected.flatten(), bins=nbins,
range=(-100, vmax), histtype='step', log=True, label = 'All') range=(-100, vmax), histtype='step', log=True, label = 'All')
_ = ax.hist(corrected[gains == 0].flatten(), bins=nbins, range=(-100, vmax), _ = ax.hist(corrected[gains == 0].flatten(), bins=nbins, range=(-100, vmax),
alpha=0.5, log=True, label='High gain', color='green') alpha=0.5, log=True, label='High gain', color='green')
_ = ax.hist(corrected[gains == 1].flatten(), bins=nbins, range=(-100, vmax), _ = ax.hist(corrected[gains == 1].flatten(), bins=nbins, range=(-100, vmax),
alpha=0.5, log=True, label='Medium gain', color='red') alpha=0.5, log=True, label='Medium gain', color='red')
_ = ax.hist(corrected[gains == 2].flatten(), bins=nbins, _ = ax.hist(corrected[gains == 2].flatten(), bins=nbins,
range=(-100, vmax), alpha=0.5, log=True, label='Low gain', color='yellow') range=(-100, vmax), alpha=0.5, log=True, label='Low gain', color='yellow')
_ = ax.legend() _ = ax.legend()
_ = ax.grid() _ = ax.grid()
_ = plt.xlabel('[ADU]') _ = plt.xlabel('[ADU]')
_ = plt.ylabel('Counts') _ = plt.ylabel('Counts')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
display(Markdown('### Maximum GAIN Preview ###\n')) display(Markdown('### Maximum GAIN Preview ###\n'))
display(Markdown(f'The per pixel maximum across one train for the digitized gain')) display(Markdown(f'The per pixel maximum across one train for the digitized gain'))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
ax = geom.plot_data_fast(np.max(gains, axis=0), ax=ax, ax = geom.plot_data_fast(np.max(gains, axis=0), ax=ax,
cmap="jet", vmin=-1, vmax=3) cmap="jet", vmin=-1, vmax=3)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Bad Pixels ## ## Bad Pixels ##
The mask contains dedicated entries for all pixels and memory cells as well as all three gains stages. Each mask entry is encoded in 32 bits as: The mask contains dedicated entries for all pixels and memory cells as well as all three gains stages. Each mask entry is encoded in 32 bits as:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
table = [] table = []
for item in BadPixels: for item in BadPixels:
table.append((item.name, "{:016b}".format(item.value))) table.append((item.name, "{:016b}".format(item.value)))
md = display(Latex(tabulate.tabulate(table, tablefmt='latex', md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
headers=["Bad pixel type", "Bit mask"]))) headers=["Bad pixel type", "Bit mask"])))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
display(Markdown(f'### Single Shot Bad Pixels ### \n')) display(Markdown(f'### Single Shot Bad Pixels ### \n'))
display(Markdown(f'A single shot bad pixel map from cell {np.max(cellId[cell_id_preview])} \n')) display(Markdown(f'A single shot bad pixel map from cell {np.max(cellId[cell_id_preview])} \n'))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
ax = geom.plot_data_fast(np.log2(mask[cell_id_preview]), ax=ax, vmin=0, vmax=32, cmap="jet") ax = geom.plot_data_fast(np.log2(mask[cell_id_preview]), ax=ax, vmin=0, vmax=32, cmap="jet")
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Percentage of Bad Pixels across one train ### ### Percentage of Bad Pixels across one train ###
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
ax = geom.plot_data_fast(np.mean(mask>0, axis=0), ax = geom.plot_data_fast(np.mean(mask>0, axis=0),
vmin=0, ax=ax, vmax=1, cmap="jet") vmin=0, ax=ax, vmax=1, cmap="jet")
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Percentage of Bad Pixels across one train. Only Dark Related ### ### Percentage of Bad Pixels across one train. Only Dark Related ###
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(20, 10)) fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
cm = np.copy(mask) cm = np.copy(mask)
cm[cm > BadPixels.NO_DARK_DATA.value] = 0 cm[cm > BadPixels.NO_DARK_DATA.value] = 0
ax = geom.plot_data_fast(np.mean(cm>0, axis=0), ax = geom.plot_data_fast(np.mean(cm>0, axis=0),
vmin=0, ax=ax, vmax=1, cmap="jet") vmin=0, ax=ax, vmax=1, cmap="jet")
``` ```
......
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# AGIPD Retrieving Constants Pre-correction # # AGIPD Retrieving Constants Pre-correction #
Author: European XFEL Detector Group, Version: 1.0 Author: European XFEL Detector Group, Version: 1.0
Retrieving Required Constants for Offline Calibration of the AGIPD Detector Retrieving Required Constants for Offline Calibration of the AGIPD Detector
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
cluster_profile = "noDB" cluster_profile = "noDB"
in_folder = "/gpfs/exfel/exp/SPB/202030/p900119/raw" # the folder to read data from, required in_folder = "/gpfs/exfel/exp/SPB/202030/p900119/raw" # the folder to read data from, required
out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/AGIPD_" # the folder to output to, required out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/AGIPD_" # the folder to output to, required
sequences = [-1] # sequences to correct, set to -1 for all, range allowed sequences = [-1] # sequences to correct, set to -1 for all, range allowed
modules = [-1] # modules to correct, set to -1 for all, range allowed modules = [-1] # modules to correct, set to -1 for all, range allowed
run = 80 # runs to process, required run = 80 # runs to process, required
karabo_id = "SPB_DET_AGIPD1M-1" # karabo karabo_id karabo_id = "SPB_DET_AGIPD1M-1" # karabo karabo_id
karabo_da = ['-1'] # a list of data aggregators names, Default [-1] for selecting all data aggregators karabo_da = ['-1'] # a list of data aggregators names, Default [-1] for selecting all data aggregators
path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data
h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP_TEST' # path to control information h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP_TEST' # path to control information
karabo_id_control = "SPB_IRU_AGIPD1M1" # karabo-id for control device karabo_id_control = "SPB_IRU_AGIPD1M1" # karabo-id for control device
karabo_da_control = 'AGIPD1MCTRL00' # karabo DA for control infromation karabo_da_control = 'AGIPD1MCTRL00' # karabo DA for control infromation
use_dir_creation_date = True # use the creation data of the input dir for database queries use_dir_creation_date = True # use the creation data of the input dir for database queries
cal_db_interface = "tcp://max-exfl016:8015#8045" # the database interface to use cal_db_interface = "tcp://max-exfl016:8015#8045" # the database interface to use
creation_date_offset = "00:00:00" # add an offset to creation date, e.g. to get different constants creation_date_offset = "00:00:00" # add an offset to creation date, e.g. to get different constants
slopes_ff_from_files = "" # Path to locally stored SlopesFF and BadPixelsFF constants
calfile = "" # path to calibration file. Leave empty if all data should come from DB calfile = "" # path to calibration file. Leave empty if all data should come from DB
nodb = False # if set only file-based constants will be used nodb = False # if set only file-based constants will be used
mem_cells = 0 # number of memory cells used, set to 0 to automatically infer mem_cells = 0 # number of memory cells used, set to 0 to automatically infer
bias_voltage = 300 bias_voltage = 300
acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine
gain_setting = 0.1 # the gain setting, use 0.1 to try to auto-determine gain_setting = 0.1 # the gain setting, use 0.1 to try to auto-determine
photon_energy = 9.2 # photon energy in keV photon_energy = 9.2 # photon energy in keV
max_cells_db_dark = 0 # set to a value different than 0 to use this value for dark data DB queries max_cells_db_dark = 0 # set to a value different than 0 to use this value for dark data DB queries
max_cells_db = 0 # set to a value different than 0 to use this value for DB queries max_cells_db = 0 # set to a value different than 0 to use this value for DB queries
# Correction Booleans # Correction Booleans
only_offset = False # Apply only Offset correction. if False, Offset is applied by Default. if True, Offset is only applied. only_offset = False # Apply only Offset correction. if False, Offset is applied by Default. if True, Offset is only applied.
rel_gain = False # do relative gain correction based on PC data rel_gain = False # do relative gain correction based on PC data
xray_gain = True # do relative gain correction based on xray data xray_gain = True # do relative gain correction based on xray data
blc_noise = False # if set, baseline correction via noise peak location is attempted blc_noise = False # if set, baseline correction via noise peak location is attempted
blc_stripes = False # if set, baseline corrected via stripes blc_stripes = False # if set, baseline corrected via stripes
blc_hmatch = False # if set, base line correction via histogram matching is attempted blc_hmatch = False # if set, base line correction via histogram matching is attempted
match_asics = False # if set, inner ASIC borders are matched to the same signal level match_asics = False # if set, inner ASIC borders are matched to the same signal level
adjust_mg_baseline = False # adjust medium gain baseline to match highest high gain value adjust_mg_baseline = False # adjust medium gain baseline to match highest high gain value
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Fill dictionaries comprising bools and arguments for correction and data analysis # Fill dictionaries comprising bools and arguments for correction and data analysis
# Here the herarichy and dependability for correction booleans are defined # Here the herarichy and dependability for correction booleans are defined
corr_bools = {} corr_bools = {}
# offset is at the bottom of AGIPD correction pyramid. # offset is at the bottom of AGIPD correction pyramid.
corr_bools["only_offset"] = only_offset corr_bools["only_offset"] = only_offset
# Dont apply any corrections if only_offset is requested # Dont apply any corrections if only_offset is requested
if not only_offset: if not only_offset:
corr_bools["adjust_mg_baseline"] = adjust_mg_baseline corr_bools["adjust_mg_baseline"] = adjust_mg_baseline
corr_bools["rel_gain"] = rel_gain corr_bools["rel_gain"] = rel_gain
corr_bools["xray_corr"] = xray_gain corr_bools["xray_corr"] = xray_gain
corr_bools["blc_noise"] = blc_noise corr_bools["blc_noise"] = blc_noise
corr_bools["blc_hmatch"] = blc_hmatch corr_bools["blc_hmatch"] = blc_hmatch
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import sys import sys
from collections import OrderedDict from collections import OrderedDict
# make sure a cluster is running with ipcluster start --n=32, give it a while to start
import os import os
import h5py import h5py
import numpy as np import numpy as np
import matplotlib import matplotlib
matplotlib.use("agg") matplotlib.use("agg")
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from ipyparallel import Client import multiprocessing as mp
print(f"Connecting to profile {cluster_profile}")
view = Client(profile=cluster_profile)[:]
view.use_dill()
from iCalibrationDB import Constants, Conditions, Detectors from iCalibrationDB import Constants, Conditions, Detectors
from cal_tools.tools import (map_modules_from_folder, get_dir_creation_date) from cal_tools.tools import (map_modules_from_folder, get_dir_creation_date)
from cal_tools.agipdlib import get_gain_setting from cal_tools.agipdlib import get_gain_setting
from dateutil import parser from dateutil import parser
from datetime import timedelta from datetime import timedelta
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
max_cells = mem_cells max_cells = mem_cells
creation_time = None creation_time = None
if use_dir_creation_date: if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run) creation_time = get_dir_creation_date(in_folder, run)
offset = parser.parse(creation_date_offset) offset = parser.parse(creation_date_offset)
delta = timedelta(hours=offset.hour, minutes=offset.minute, seconds=offset.second) delta = timedelta(hours=offset.hour, minutes=offset.minute, seconds=offset.second)
creation_time += delta creation_time += delta
print(f"Using {creation_time} as creation time") print(f"Using {creation_time} as creation time")
if sequences[0] == -1: if sequences[0] == -1:
sequences = None sequences = None
if in_folder[-1] == "/": if in_folder[-1] == "/":
in_folder = in_folder[:-1] in_folder = in_folder[:-1]
print(f"Outputting to {out_folder}") print(f"Outputting to {out_folder}")
os.makedirs(out_folder, exist_ok=True) os.makedirs(out_folder, exist_ok=True)
import warnings import warnings
warnings.filterwarnings('ignore') warnings.filterwarnings('ignore')
from cal_tools.agipdlib import SnowResolution from cal_tools.agipdlib import SnowResolution
melt_snow = False if corr_bools["only_offset"] else SnowResolution.NONE melt_snow = False if corr_bools["only_offset"] else SnowResolution.NONE
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
control_fname = f'{in_folder}/r{run:04d}/RAW-R{run:04d}-{karabo_da_control}-S00000.h5' control_fname = f'{in_folder}/r{run:04d}/RAW-R{run:04d}-{karabo_da_control}-S00000.h5'
h5path_ctrl = h5path_ctrl.format(karabo_id_control) h5path_ctrl = h5path_ctrl.format(karabo_id_control)
if gain_setting == 0.1: if gain_setting == 0.1:
if creation_time.replace(tzinfo=None) < parser.parse('2020-01-31'): if creation_time.replace(tzinfo=None) < parser.parse('2020-01-31'):
print("Set gain-setting to None for runs taken before 2020-01-31") print("Set gain-setting to None for runs taken before 2020-01-31")
gain_setting = None gain_setting = None
else: else:
try: try:
gain_setting = get_gain_setting(control_fname, h5path_ctrl) gain_setting = get_gain_setting(control_fname, h5path_ctrl)
except Exception as e: except Exception as e:
print(f'ERROR: while reading gain setting from: \n{control_fname}') print(f'ERROR: while reading gain setting from: \n{control_fname}')
print(e) print(e)
print("Set gain setting to 0") print("Set gain setting to 0")
gain_setting = 0 gain_setting = 0
print(f"Gain setting: {gain_setting}") print(f"Gain setting: {gain_setting}")
print(f"Detector in use is {karabo_id}") print(f"Detector in use is {karabo_id}")
# Extracting Instrument string # Extracting Instrument string
instrument = karabo_id.split("_")[0] instrument = karabo_id.split("_")[0]
# Evaluate detector instance for mapping # Evaluate detector instance for mapping
if instrument == "SPB": if instrument == "SPB":
dinstance = "AGIPD1M1" dinstance = "AGIPD1M1"
nmods = 16 nmods = 16
elif instrument == "MID": elif instrument == "MID":
dinstance = "AGIPD1M2" dinstance = "AGIPD1M2"
nmods = 16 nmods = 16
# TODO: Remove DETLAB # TODO: Remove DETLAB
elif instrument == "HED" or instrument == "DETLAB": elif instrument == "HED" or instrument == "DETLAB":
dinstance = "AGIPD500K" dinstance = "AGIPD500K"
nmods = 8 nmods = 8
print(f"Instrument {instrument}") print(f"Instrument {instrument}")
print(f"Detector instance {dinstance}") print(f"Detector instance {dinstance}")
if karabo_da[0] == '-1': if karabo_da[0] == '-1':
if modules[0] == -1: if modules[0] == -1:
modules = list(range(nmods)) modules = list(range(nmods))
karabo_da = ["AGIPD{:02d}".format(i) for i in modules] karabo_da = ["AGIPD{:02d}".format(i) for i in modules]
else: else:
modules = [int(x[-2:]) for x in karabo_da] modules = [int(x[-2:]) for x in karabo_da]
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# set everything up filewise # set everything up filewise
print(f"Checking the files before retrieving constants") print(f"Checking the files before retrieving constants")
mmf = map_modules_from_folder(in_folder, run, path_template, karabo_da, sequences) mmf = map_modules_from_folder(in_folder, run, path_template, karabo_da, sequences)
mapped_files, mod_ids, total_sequences, sequences_qm, _ = mmf mapped_files, mod_ids, total_sequences, sequences_qm, _ = mmf
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Retrieve Constants ## ## Retrieve Constants ##
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from functools import partial from functools import partial
import yaml import yaml
def retrieve_constants(karabo_id, bias_voltage, max_cells, acq_rate, def retrieve_constants(karabo_id, bias_voltage, max_cells, acq_rate,
gain_setting, photon_energy, only_dark, nodb_with_dark, gain_setting, photon_energy, only_dark, nodb_with_dark,
cal_db_interface, creation_time, cal_db_interface, creation_time,
corr_bools, pc_bools, inp): corr_bools, pc_bools, inp):
""" """
Retreive constant for each module in parallel and produce a dictionary Retreive constant for each module in parallel and produce a dictionary
with the creation-time and constant file path. with the creation-time and constant file path.
:param karabo_id: (STR) Karabo ID :param karabo_id: (STR) Karabo ID
:param bias_voltage: (FLOAT) Bias Voltage :param bias_voltage: (FLOAT) Bias Voltage
:param max_cells: (INT) Memory cells :param max_cells: (INT) Memory cells
:param acq_rate: (FLOAT) Acquisition Rate :param acq_rate: (FLOAT) Acquisition Rate
:param gain_setting: (FLOAT) Gain setting :param gain_setting: (FLOAT) Gain setting
:param photon_energy: (FLOAT) Photon Energy :param photon_energy: (FLOAT) Photon Energy
:param only_dark: (BOOL) only retrieve dark constants :param only_dark: (BOOL) only retrieve dark constants
:param nodb_with_dark: (BOOL) no constant retrieval even for dark :param nodb_with_dark: (BOOL) no constant retrieval even for dark
:param cal_db_interface: (STR) the database interface port :param cal_db_interface: (STR) the database interface port
:param creation_time: (STR) raw data creation time :param creation_time: (STR) raw data creation time
:param corr_bools: (DICT) A dictionary with bools for applying requested corrections :param corr_bools: (DICT) A dictionary with bools for applying requested corrections
:param pc_bools: (LIST) list of bools to retrieve pulse capacitor constants :param pc_bools: (LIST) list of bools to retrieve pulse capacitor constants
:param inp: (LIST) input for the parallel cluster of the partial function :param inp: (LIST) input for the parallel cluster of the partial function
:return: :return:
mdata_dict: (DICT) dictionary with the metadata for the retrieved constants mdata_dict: (DICT) dictionary with the metadata for the retrieved constants
dev.device_name: (STR) device name dev.device_name: (STR) device name
""" """
import numpy as np import numpy as np
import sys import sys
import traceback import traceback
from cal_tools.agipdlib import get_num_cells, get_acq_rate from cal_tools.agipdlib import get_num_cells, get_acq_rate
from cal_tools.agipdutils import assemble_constant_dict from cal_tools.agipdutils import assemble_constant_dict
from cal_tools.tools import get_from_db from cal_tools.tools import get_from_db
from iCalibrationDB import Constants, Conditions, Detectors from iCalibrationDB import Constants, Conditions, Detectors
err = None err = None
qm_files, qm, dev, idx = inp qm_files, qm, dev, idx = inp
# get number of memory cells from a sequence file with image data # get number of memory cells from a sequence file with image data
for f in qm_files: for f in qm_files:
if not max_cells: if not max_cells:
max_cells = get_num_cells(f, karabo_id, idx) max_cells = get_num_cells(f, karabo_id, idx)
if max_cells is None: if max_cells is None:
if f != qm_files[-1]: if f != qm_files[-1]:
continue continue
else: else:
raise ValueError(f"No raw images found for {qm} for all sequences") raise ValueError(f"No raw images found for {qm} for all sequences")
else: else:
cells = np.arange(max_cells) cells = np.arange(max_cells)
# get out of the loop, # get out of the loop,
# if max_cells is successfully calculated. # if max_cells is successfully calculated.
break break
if acq_rate == 0.: if acq_rate == 0.:
acq_rate = get_acq_rate((f, karabo_id, idx)) acq_rate = get_acq_rate((f, karabo_id, idx))
# avoid creating retireving constant, if requested. # avoid creating retireving constant, if requested.
if not nodb_with_dark: if not nodb_with_dark:
const_dict = assemble_constant_dict(corr_bools, pc_bools, max_cells, bias_voltage, const_dict = assemble_constant_dict(corr_bools, pc_bools, max_cells, bias_voltage,
gain_setting, acq_rate, photon_energy, gain_setting, acq_rate, photon_energy,
beam_energy=None, only_dark=only_dark) beam_energy=None, only_dark=only_dark)
# Retrieve multiple constants through an input dictionary # Retrieve multiple constants through an input dictionary
# to return a dict of useful metadata. # to return a dict of useful metadata.
mdata_dict = dict() mdata_dict = dict()
for cname, cval in const_dict.items(): for cname, cval in const_dict.items():
try: # saving metadata in a dict
condition = getattr(Conditions, cval[2][0]).AGIPD(**cval[2][1]) mdata_dict[cname] = dict()
co, mdata = \
get_from_db(dev, getattr(Constants.AGIPD, cname)(), if slopes_ff_from_files and cname in ["SlopesFF", "BadPixelsFF"]:
condition, getattr(np, cval[0])(cval[1]), mdata_dict[cname]["file-path"] = f"{slopes_ff_from_files}/slopesff_bpmask_module_{qm}.h5"
cal_db_interface, creation_time, meta_only=True, verbosity=0) mdata_dict[cname]["creation-time"] = "00:00:00"
mdata_const = mdata.calibration_constant_version else:
# saving metadata in a dict try:
mdata_dict[cname] = dict() condition = getattr(Conditions, cval[2][0]).AGIPD(**cval[2][1])
# check if constant was sucessfully retrieved. co, mdata = \
if mdata.comm_db_success: get_from_db(dev, getattr(Constants.AGIPD, cname)(),
mdata_dict[cname]["file-path"] = f"{mdata_const.hdf5path}" \ condition, getattr(np, cval[0])(cval[1]),
f"{mdata_const.filename}" cal_db_interface, creation_time, meta_only=True, verbosity=0)
mdata_dict[cname]["creation-time"] = f"{mdata_const.begin_at}" mdata_const = mdata.calibration_constant_version
else:
mdata_dict[cname]["file-path"] = const_dict[cname][:2] # check if constant was sucessfully retrieved.
mdata_dict[cname]["creation-time"] = None if mdata.comm_db_success:
except Exception as e: mdata_dict[cname]["file-path"] = f"{mdata_const.hdf5path}" \
err = f"Error: {e}, Traceback: {traceback.format_exc()}" f"{mdata_const.filename}"
print(err) mdata_dict[cname]["creation-time"] = f"{mdata_const.begin_at}"
else:
mdata_dict[cname]["file-path"] = const_dict[cname][:2]
mdata_dict[cname]["creation-time"] = None
except Exception as e:
err = f"Error: {e}, Traceback: {traceback.format_exc()}"
print(err)
return qm, mdata_dict, dev.device_name, acq_rate, max_cells, err return qm, mdata_dict, dev.device_name, acq_rate, max_cells, err
pc_bools = [corr_bools.get("rel_gain"), pc_bools = [corr_bools.get("rel_gain"),
corr_bools.get("adjust_mg_baseline"), corr_bools.get("adjust_mg_baseline"),
corr_bools.get('blc_noise'), corr_bools.get('blc_noise'),
corr_bools.get('blc_hmatch'), corr_bools.get('blc_hmatch'),
corr_bools.get('blc_stripes'), corr_bools.get('blc_stripes'),
melt_snow] melt_snow]
inp = [] inp = []
only_dark = False only_dark = False
nodb_with_dark = False nodb_with_dark = False
if not nodb: if not nodb:
only_dark=(calfile != "") only_dark=(calfile != "")
if calfile != "" and not corr_bools["only_offset"]: if calfile != "" and not corr_bools["only_offset"]:
nodb_with_dark = nodb nodb_with_dark = nodb
# A dict to connect virtual device # A dict to connect virtual device
# to actual device name. # to actual device name.
for i in modules: for i in modules:
qm = f"Q{i//4+1}M{i%4+1}" qm = f"Q{i//4+1}M{i%4+1}"
if qm in mapped_files and not mapped_files[qm].empty(): if qm in mapped_files and not mapped_files[qm].empty():
device = getattr(getattr(Detectors, dinstance), qm) device = getattr(getattr(Detectors, dinstance), qm)
qm_files = [str(mapped_files[qm].get()) for _ in range(mapped_files[qm].qsize())] qm_files = [str(mapped_files[qm].get()) for _ in range(mapped_files[qm].qsize())]
else: else:
print(f"Skipping {qm}") print(f"Skipping {qm}")
continue continue
inp.append((qm_files, qm, device, i)) inp.append((qm_files, qm, device, i))
p = partial(retrieve_constants, karabo_id, bias_voltage, max_cells, p = partial(retrieve_constants, karabo_id, bias_voltage, max_cells,
acq_rate, gain_setting, photon_energy, only_dark, nodb_with_dark, acq_rate, gain_setting, photon_energy, only_dark, nodb_with_dark,
cal_db_interface, creation_time, cal_db_interface, creation_time,
corr_bools, pc_bools) corr_bools, pc_bools)
results = view.map_sync(p, inp) with mp.Pool(processes=16) as pool:
#results = list(map(p, inp)) results = pool.map(p, inp)
mod_dev = dict() mod_dev = dict()
mdata_dict = dict() mdata_dict = dict()
for r in results: for r in results:
if r: if r:
qm, md_dict, dname, acq_rate, max_cells, err = r qm, md_dict, dname, acq_rate, max_cells, err = r
mod_dev[dname] = {"mod": qm, "err": err} mod_dev[dname] = {"mod": qm, "err": err}
if err: if err:
print(f"Error for module {qm}: {err}") print(f"Error for module {qm}: {err}")
mdata_dict[dname] = md_dict mdata_dict[dname] = md_dict
# check if it is requested not to retrieve any constants from the database # check if it is requested not to retrieve any constants from the database
if not nodb_with_dark: if not nodb_with_dark:
with open(f"{out_folder}/retrieved_constants.yml", "w") as outfile: with open(f"{out_folder}/retrieved_constants.yml", "w") as outfile:
yaml.safe_dump(mdata_dict, outfile) yaml.safe_dump(mdata_dict, outfile)
print("\nRetrieved constants for modules: ", print("\nRetrieved constants for modules: ",
f"{[', '.join([f'Q{x//4+1}M{x%4+1}' for x in modules])]}") f"{[', '.join([f'Q{x//4+1}M{x%4+1}' for x in modules])]}")
print(f"Operating conditions are:\n• Bias voltage: {bias_voltage}\n• Memory cells: {max_cells}\n" print(f"Operating conditions are:\n• Bias voltage: {bias_voltage}\n• Memory cells: {max_cells}\n"
f"• Acquisition rate: {acq_rate}\n• Gain setting: {gain_setting}\n• Photon Energy: {photon_energy}\n") f"• Acquisition rate: {acq_rate}\n• Gain setting: {gain_setting}\n• Photon Energy: {photon_energy}\n")
print(f"Constant metadata is saved in retrieved_constants.yml\n") print(f"Constant metadata is saved in retrieved_constants.yml\n")
else: else:
print("No constants were retrieved as calibrated files will be used.") print("No constants were retrieved as calibrated files will be used.")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
print("Constants are retrieved with creation time: ") print("Constants are retrieved with creation time: ")
i = 0 i = 0
when = dict() when = dict()
to_store = [] to_store = []
for dname, dinfo in mod_dev.items(): for dname, dinfo in mod_dev.items():
print(dinfo["mod"], ":") print(dinfo["mod"], ":")
line = [dinfo["mod"]] line = [dinfo["mod"]]
if dname in mdata_dict: if dname in mdata_dict:
for cname, mdata in mdata_dict[dname].items(): for cname, mdata in mdata_dict[dname].items():
if hasattr(mdata["creation-time"], 'strftime'): if hasattr(mdata["creation-time"], 'strftime'):
mdata["creation-time"] = mdata["creation-time"].strftime('%y-%m-%d %H:%M') mdata["creation-time"] = mdata["creation-time"].strftime('%y-%m-%d %H:%M')
print(f'{cname:.<12s}', mdata["creation-time"]) print(f'{cname:.<12s}', mdata["creation-time"])
# Store few time stamps if exists # Store few time stamps if exists
# Add NA to keep array structure # Add NA to keep array structure
for cname in ['Offset', 'SlopesPC', 'SlopesFF']: for cname in ['Offset', 'SlopesPC', 'SlopesFF']:
if not dname in mdata_dict or dinfo["err"]: if not dname in mdata_dict or dinfo["err"]:
line.append('Err') line.append('Err')
else: else:
if cname in mdata_dict[dname]: if cname in mdata_dict[dname]:
if mdata_dict[dname][cname]["creation-time"]: if mdata_dict[dname][cname]["creation-time"]:
line.append(mdata_dict[dname][cname]["creation-time"]) line.append(mdata_dict[dname][cname]["creation-time"])
else: else:
line.append('NA') line.append('NA')
else: else:
line.append('NA') line.append('NA')
to_store.append(line) to_store.append(line)
i += 1 i += 1
if sequences: if sequences:
seq_num = sequences[0] seq_num = sequences[0]
else: else:
# if sequences[0] changed to None as it was -1 # if sequences[0] changed to None as it was -1
seq_num = 0 seq_num = 0
with open(f"{out_folder}/retrieved_constants.yml","r") as fyml: with open(f"{out_folder}/retrieved_constants.yml","r") as fyml:
time_summary = yaml.safe_load(fyml) time_summary = yaml.safe_load(fyml)
time_summary.update({"time-summary": { time_summary.update({"time-summary": {
"SAll":to_store "SAll":to_store
}}) }})
with open(f"{out_folder}/retrieved_constants.yml","w") as fyml: with open(f"{out_folder}/retrieved_constants.yml","w") as fyml:
yaml.safe_dump(time_summary, fyml) yaml.safe_dump(time_summary, fyml)
``` ```
%% Cell type:code id: tags:
``` python
```
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
import numpy as np
from cal_tools.agipdutils_ff import get_mask, set_par_limits
def test_get_mask():
fit_summary = {
'chi2_ndof': 1.674524751845516,
'g0n': 6031.641198873036,
'error_g0n': 94.63055028459667,
'limit_g0n': np.array([0.0, None]),
'fix_g0n': False,
'g0mean': -13.711814669099589,
'error_g0mean': 0.2532017427306297,
'limit_g0mean': np.array([-30, 30]),
'fix_g0mean': False,
'g0sigma': 13.478502058651568,
'error_g0sigma': 0.2588135637661919,
'limit_g0sigma': np.array([0, 30]),
'fix_g0sigma': False,
'g1n': 4337.126861254491,
'error_g1n': 69.764180118274,
'limit_g1n': np.array([0, None]),
'fix_g1n': False,
'g1mean': 53.90265411499657,
'error_g1mean': 0.27585684670864746,
'limit_g1mean': None,
'fix_g1mean': False,
'g1sigma': 15.687448834904817,
'error_g1sigma': 0.2951166525483524,
'limit_g1sigma': np.array([0, 35]),
'fix_g1sigma': False,
'g2n': 1542.531700635003,
'error_g2n': 43.20145030604567,
'limit_g2n': np.array([0, None]),
'fix_g2n': False,
'g2mean': 120.98535387591575,
'error_g2mean': 0.509566354942716,
'limit_g2mean': None,
'fix_g2mean': False,
'g2sigma': 15.550452880533143,
'error_g2sigma': 0.5003254358001863,
'limit_g2sigma': np.array([0, 40]),
'fix_g2sigma': False,
'g3n': 1261189.2282171287,
'error_g3n': 1261190.2282163086,
'limit_g3n': np.array([0, None]),
'fix_g3n': False,
'g3mean': 348.68766379647343,
'error_g3mean': 17.23872285713375,
'limit_g3mean': None,
'fix_g3mean': False,
'g3sigma': 44.83987230934497,
'error_g3sigma': 30.956164693249242,
'limit_g3sigma': np.array([0, 45]),
'fix_g3sigma': False,
'fval': 336.5794751209487,
'edm': 0.00011660826330754263,
'tolerance': 0.1,
'nfcn': 4620,
'ncalls': 4620,
'up': 1.0,
'is_valid': True,
'has_valid_parameters': True,
'has_accurate_covar': True,
'has_posdef_covar': True,
'has_made_posdef_covar': False,
'hesse_failed': False,
'has_covariance': True,
'is_above_max_edm': False,
'has_reached_call_limit': False}
peak_lim = [-30, 30]
d0_lim = [10, 80]
chi2_lim = [0, 3.0]
peak_width_lim = np.array([[0.9, 1.55], [0.95, 1.65]])
mask = get_mask(fit_summary, peak_lim, d0_lim, chi2_lim, peak_width_lim)
assert mask == 0
def test_set_par_limits():
peak_range = np.array([[-30, 30],
[35, 70],
[95, 135],
[145, 220]])
peak_norm_range = np.array([[0.0, None],
[0, None],
[0, None],
[0, None]])
peak_width_range = np.array([[0, 30],
[0, 35],
[0, 40],
[0, 45]])
parameters = {
'g0sigma': 9.620186459204016,
'g0n': 5659.0,
'g0mean': -3.224686340342817,
'g1sigma': 8.149415371586683,
'g1n': 3612.0,
'g1mean': 54.6281838316722,
'g2sigma': 9.830124777667839,
'g2n': 1442.0,
'g2mean': 114.92510402219139,
'g3sigma': 15.336595220228498,
'g3n': 474.0,
'g3mean': 167.0295358649789}
expected = {
'g0sigma': 9.620186459204016,
'g0n': 5659.0,
'g0mean': -3.224686340342817,
'g1sigma': 8.149415371586683,
'g1n': 3612.0,
'g1mean': 54.6281838316722,
'g2sigma': 9.830124777667839,
'g2n': 1442.0,
'g2mean': 114.92510402219139,
'g3sigma': 15.336595220228498,
'g3n': 474.0,
'g3mean': 167.0295358649789,
'limit_g0n': np.array([0.0, None]),
'limit_g0mean': np.array([-30, 30]),
'limit_g0sigma': np.array([0, 30]),
'limit_g1n': np.array([0, None]),
'limit_g1mean': np.array([35, 70]),
'limit_g1sigma': np.array([0, 35]),
'limit_g2n': np.array([0, None]),
'limit_g2mean': np.array([95, 135]),
'limit_g2sigma': np.array([0, 40]),
'limit_g3n': np.array([0, None]),
'limit_g3mean': np.array([145, 220]),
'limit_g3sigma': np.array([0, 45])}
set_par_limits(parameters, peak_range, peak_norm_range, peak_width_range)
assert parameters.keys() == expected.keys()
for key in parameters.keys():
if isinstance(parameters[key], np.ndarray):
assert np.all(parameters[key] == expected[key])
else:
assert parameters[key] == expected[key]
...@@ -21,6 +21,8 @@ notebooks = { ...@@ -21,6 +21,8 @@ notebooks = {
"FF": { "FF": {
"notebook": "notebook":
"notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb", "notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb",
"dep_notebooks": [
"notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_Summary.ipynb"],
"concurrency": {"parameter": "modules", "concurrency": {"parameter": "modules",
"default concurrency": 16, "default concurrency": 16,
"cluster cores": 16}, "cluster cores": 16},
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment