diff --git a/notebooks/LPD/LPD_Correct_Fast.ipynb b/notebooks/LPD/LPD_Correct_Fast.ipynb
index 04d35400aafeaa9197d79ff6a11439d697c84df0..539ce31b7c14bd5657275d2856bb99d5078a8a21 100644
--- a/notebooks/LPD/LPD_Correct_Fast.ipynb
+++ b/notebooks/LPD/LPD_Correct_Fast.ipynb
@@ -33,6 +33,7 @@
"karabo_id = 'FXE_DET_LPD1M-1' # Karabo domain for detector.\n",
"input_source = '{karabo_id}/DET/{module_index}CH0:xtdf' # Input fast data source.\n",
"output_source = '{karabo_id}/CORR/{module_index}CH0:output' # Output fast data source, empty to use same as input.\n",
+ "control_source = '{karabo_id}/COMP/FEM_MDL_COMP' # Control data source.\n",
"xgm_source = 'SA1_XTD2_XGM/DOOCS/MAIN'\n",
"xgm_pulse_count_key = 'pulseEnergy.numberOfSa1BunchesActual'\n",
"\n",
@@ -55,6 +56,9 @@
"rel_gain = True # Gain correction based on RelativeGain constant.\n",
"ff_map = True # Gain correction based on FFMap constant.\n",
"gain_amp_map = True # Gain correction based on GainAmpMap constant.\n",
+ "combine_parallel_gain = True # Combine parallel gain images into a single frame.\n",
+ "threshold_sigma_high = 5.0 # Sigma level for threshold between high and medium gain.\n",
+ "threshold_sigma_mid = 100.0 # Sigma level for threshold between medium and low gain.\n",
"\n",
"# Output options\n",
"ignore_no_frames_no_pulses = False # Whether to run without SA1 pulses AND frames.\n",
@@ -244,13 +248,24 @@
"source": [
"start = perf_counter()\n",
"\n",
+ "raw_data = xd.RunDirectory(run_folder)\n",
+ "\n",
+ "try:\n",
+ " parallel_gain = bool(raw_data[control_source.format(karabo_id=karabo_id)].run_value('femAsicGainOverride'))\n",
+ "except KeyError:\n",
+ " warning('Missing femAsicGainOverride property FEM control device, assuming auto gain')\n",
+ " parallel_gain = False\n",
+ "print('Parallel gain mode:', parallel_gain)\n",
+ "\n",
"cell_ids_pattern_s = None\n",
"if use_cell_order != 'never':\n",
+ " mem_cell_pattern = get_mem_cell_pattern(raw_data, det_inp_sources)\n",
+ " \n",
+ " if parallel_gain:\n",
+ " mem_cell_pattern = mem_cell_pattern[:len(mem_cell_pattern) // 3]\n",
+ " \n",
" # Read the order of memory cells used\n",
- " raw_data = xd.DataCollection.from_paths([e[1] for e in data_to_process])\n",
- " cell_ids_pattern_s = make_cell_order_condition(\n",
- " use_cell_order, get_mem_cell_pattern(raw_data, det_inp_sources)\n",
- " )\n",
+ " cell_ids_pattern_s = make_cell_order_condition(use_cell_order, mem_cell_pattern)\n",
"print(\"Memory cells order:\", cell_ids_pattern_s)\n",
"\n",
"conditions = LPDConditions(\n",
@@ -259,6 +274,7 @@
" feedback_capacitor=capacitor,\n",
" source_energy=photon_energy,\n",
" memory_cell_order=cell_ids_pattern_s,\n",
+ " parallel_gain=parallel_gain,\n",
" category=category,\n",
")\n",
"\n",
@@ -269,6 +285,8 @@
" expected_constants.update(['FFMap', 'BadPixelsFF'])\n",
"if gain_amp_map:\n",
" expected_constants.add('GainAmpMap')\n",
+ "if parallel_gain:\n",
+ " expected_constants.add('Noise')\n",
"\n",
"lpd_consts = CalibrationData.from_condition(\n",
" conditions,\n",
@@ -345,6 +363,7 @@
"source": [
"# These are intended in order cell, X, Y, gain\n",
"ccv_offsets = {}\n",
+ "ccv_noise = {}\n",
"ccv_gains = {}\n",
"ccv_masks = {}\n",
"\n",
@@ -352,6 +371,7 @@
"\n",
"constant_order = {\n",
" 'Offset': (2, 1, 0, 3),\n",
+ " 'Noise': (2, 1, 0, 3),\n",
" 'BadPixelsDark': (2, 1, 0, 3),\n",
" 'RelativeGain': (2, 0, 1, 3),\n",
" 'FFMap': (2, 0, 1, 3),\n",
@@ -378,6 +398,11 @@
" ccv_offsets[aggregator] = np.zeros(ccv_shape, dtype=np.float32)\n",
" \n",
" ccv_gains[aggregator] = np.ones(ccv_shape, dtype=np.float32)\n",
+ "\n",
+ " if parallel_gain and 'Noise' in consts:\n",
+ " ccv_noise[aggregator] = _prepare_data('Noise', np.float32)\n",
+ " else:\n",
+ " ccv_noise[aggregator] = np.zeros(ccv_shape, dtype=np.float32)\n",
" \n",
" if 'BadPixelsDark' in consts:\n",
" ccv_masks[aggregator] = _prepare_data('BadPixelsDark', np.uint32)\n",
@@ -433,24 +458,67 @@
" in_raw = inp_source['image.data'].ndarray().reshape(-1, 256, 256)\n",
" in_cell = inp_source['image.cellId'].ndarray().reshape(-1)\n",
" in_pulse = inp_source['image.pulseId'].ndarray().reshape(-1)\n",
+ " frame_counts = inp_source['image.data'].data_counts(labelled=False).astype(np.int32)\n",
" read_time = perf_counter() - start\n",
+ "\n",
+ " parallel_gain_indices = None\n",
+ " \n",
+ " if parallel_gain:\n",
+ " actual_frame_counts = frame_counts // 3\n",
+ "\n",
+ " # Indices map where to find each of the high/medium/low gain images for each actual\n",
+ " # frame event.\n",
+ " parallel_gain_indices = np.zeros((actual_frame_counts.sum(), 3), dtype=np.int32)\n",
+ "\n",
+ " # Build indices for high gain as a range in each train, running from the cumulative sum\n",
+ " # of apparent frames from all trains before to the actual number of frames in this train.\n",
+ " np.concatenate(\n",
+ " [np.arange(actual_frame_counts[0])] +\n",
+ " [np.arange(offset, offset+actual_count) for offset, actual_count\n",
+ " in zip(np.cumsum(frame_counts), actual_frame_counts[1:])],\n",
+ " out=parallel_gain_indices[:, 0])\n",
+ "\n",
+ " # The delta between the gain stages is the number of actual frames.\n",
+ " gain_index_deltas = np.repeat(actual_frame_counts, actual_frame_counts)\n",
+ "\n",
+ " # Build indices for medium gain and high gain by adding the gain index deltas in between\n",
+ " # each of them.\n",
+ " np.add(parallel_gain_indices[:, 0], gain_index_deltas, out=parallel_gain_indices[:, 1])\n",
+ " np.add(parallel_gain_indices[:, 1], gain_index_deltas, out=parallel_gain_indices[:, 2])\n",
+ "\n",
+ " # Pick cell and pulse IDs from high gain. This is also done if frames are not combined\n",
+ " # in order to correct corrupt tables in medium and low gain, and if needed brought back\n",
+ " # to the original shape further below.\n",
+ " in_cell = np.take(in_cell, parallel_gain_indices[:, 0])\n",
+ " in_pulse = np.take(in_pulse, parallel_gain_indices[:, 0])\n",
+ "\n",
+ " if combine_parallel_gain:\n",
+ " # Replace supposed frame counts by actual frame counts.\n",
+ " frame_counts = actual_frame_counts\n",
+ " else:\n",
+ " # Replicate corrected cell and pulse IDs from high gain to other gains.\n",
+ " in_cell = np.repeat(in_cell, actual_frame_counts)\n",
+ " in_pulse = np.repeat(in_pulse, actual_frame_counts)\n",
+ " \n",
+ " # Disable gain indices to not combine.\n",
+ " parallel_gain_indices = None\n",
" \n",
" # Allocate output arrays.\n",
- " out_data = np.zeros((in_raw.shape[0], 256, 256), dtype=np.float32)\n",
- " out_gain = np.zeros((in_raw.shape[0], 256, 256), dtype=np.uint8)\n",
- " out_mask = np.zeros((in_raw.shape[0], 256, 256), dtype=np.uint32)\n",
+ " num_frames = frame_counts.sum()\n",
+ " out_data = np.zeros((num_frames, 256, 256), dtype=np.float32)\n",
+ " out_gain = np.zeros((num_frames, 256, 256), dtype=np.uint8)\n",
+ " out_mask = np.zeros((num_frames, 256, 256), dtype=np.uint32)\n",
" \n",
" start = perf_counter()\n",
" correct_lpd_frames(in_raw, in_cell,\n",
" out_data, out_gain, out_mask,\n",
- " ccv_offsets[aggregator], ccv_gains[aggregator], ccv_masks[aggregator],\n",
- " num_threads=num_threads_per_worker)\n",
+ " ccv_offsets[aggregator], ccv_noise[aggregator], ccv_gains[aggregator], ccv_masks[aggregator],\n",
+ " parallel_gain_indices, threshold_sigma_high, threshold_sigma_mid,\n",
+ " num_threads=16)\n",
" correct_time = perf_counter() - start\n",
" \n",
- " image_counts = inp_source['image.data'].data_counts(labelled=False)\n",
- " \n",
" start = perf_counter()\n",
- " if (not outp_path.exists() or overwrite) and image_counts.sum() > 0:\n",
+ " if (not outp_path.exists() or overwrite) and num_frames > 0:\n",
" outp_source_name = output_source.format(karabo_id=karabo_id, module_index=module_index)\n",
"\n",
" with DataFile(outp_path, 'w') as outp_file: \n",
@@ -461,7 +529,7 @@
" \n",
" outp_source = outp_file.create_instrument_source(outp_source_name)\n",
" \n",
- " outp_source.create_index(image=image_counts)\n",
+ " outp_source.create_index(image=frame_counts)\n",
" outp_source.create_key('image.cellId', data=in_cell,\n",
" chunks=(min(chunks_ids, in_cell.shape[0]),))\n",
" outp_source.create_key('image.pulseId', data=in_pulse,\n",
@@ -477,11 +545,11 @@
" write_time = perf_counter() - start\n",
" \n",
" total_time = open_time + read_time + correct_time + write_time\n",
- " frame_rate = in_raw.shape[0] / total_time\n",
+ " frame_rate = num_frames / total_time\n",
" \n",
" print('{}\\t{}\\t{:.3f}\\t{:.3f}\\t{:.3f}\\t{:.3f}\\t{:.3f}\\t{}\\t{:.1f}'.format(\n",
" wid, aggregator, open_time, read_time, correct_time, write_time, total_time,\n",
- " in_raw.shape[0], frame_rate))\n",
+ " num_frames, frame_rate))\n",
" \n",
" in_raw = None\n",
" in_cell = None\n",
diff --git a/src/cal_tools/lpdalgs.pyx b/src/cal_tools/lpdalgs.pyx
index 66b8b097c71a4a7cb6b1681e4678e0f04c2f43e7..6c746172de4cfc109d20c0bf3e3c780ee35e9d3d 100644
--- a/src/cal_tools/lpdalgs.pyx
+++ b/src/cal_tools/lpdalgs.pyx
@@ -1,13 +1,14 @@
+
from cython cimport boundscheck, wraparound, cdivision
from cython.view cimport contiguous
from cython.parallel cimport prange
from cal_tools.enums import BadPixels
-ctypedef unsigned short cell_t
ctypedef unsigned short raw_t
ctypedef float data_t
ctypedef unsigned char gain_t
+ctypedef unsigned short cell_t
ctypedef unsigned int mask_t
cdef mask_t WRONG_GAIN_VALUE = BadPixels.WRONG_GAIN_VALUE, \
@@ -17,7 +18,6 @@ cdef mask_t WRONG_GAIN_VALUE = BadPixels.WRONG_GAIN_VALUE, \
@boundscheck(False)
@wraparound(False)
-@cdivision(True)
def correct_lpd_frames(
# (frame, x, y)
raw_t[:, :, ::contiguous] in_raw,
@@ -29,10 +29,17 @@ def correct_lpd_frames(
mask_t[:, :, ::contiguous] out_mask,
# (cell, x, y, gain)
- float[:, :, :, ::contiguous] ccv_offset,
- float[:, :, :, ::contiguous] ccv_gain,
+ data_t[:, :, :, ::contiguous] ccv_offset,
+ data_t[:, :, :, ::contiguous] ccv_noise,
+ data_t[:, :, :, ::contiguous] ccv_gain,
mask_t[:, :, :, ::contiguous] ccv_mask,
+ # (frame, gain)
+ int[:, ::contiguous] parallel_gain_indices,
+
+ data_t threshold_sigma_high,
+ data_t threshold_sigma_mid,
+
int num_threads=1,
):
cdef int frame, ss, fs
@@ -41,25 +48,34 @@ def correct_lpd_frames(
cdef gain_t gain
cdef mask_t mask
- for frame in prange(in_raw.shape[0], nogil=True, num_threads=num_threads):
+ cdef bint adaptive_gain = parallel_gain_indices is None
+
+ for frame in prange(out_data.shape[0], nogil=True, num_threads=num_threads):
cell = in_cell[frame]
for ss in range(in_raw.shape[1]):
for fs in range(in_raw.shape[2]):
- # Decode intensity and gain from raw data.
- data = <data_t>(in_raw[frame, ss, fs] & 0xFFF)
- gain = <gain_t>((in_raw[frame, ss, fs] & 0x3000) >> 12)
+ if adaptive_gain:
+ # Decode intensity and gain from raw data.
+ data = <data_t>(in_raw[frame, ss, fs] & 0xFFF)
+ gain = <gain_t>((in_raw[frame, ss, fs] & 0x3000) >> 12)
+
+ else:
+ _parallel_gain_thresholding(
+ in_raw, parallel_gain_indices, ccv_noise,
+ threshold_sigma_high, threshold_sigma_mid,
+ frame, cell, ss, fs,
+ &data, &gain)
if gain <= 2:
+ data = data - ccv_offset[cell, ss, fs, gain]
+ data = data * ccv_gain[cell, ss, fs, gain]
mask = ccv_mask[cell, ss, fs, gain]
else:
data = 0.0
gain = 0
mask = WRONG_GAIN_VALUE
- data = data - ccv_offset[cell, ss, fs, gain]
- data = data * ccv_gain[cell, ss, fs, gain]
-
if data > 1e7 or data < -1e7:
data = 0.0
mask = mask | VALUE_OUT_OF_RANGE
@@ -67,3 +83,58 @@ def correct_lpd_frames(
out_data[frame, ss, fs] = data
out_gain[frame, ss, fs] = gain
out_mask[frame, ss, fs] = mask
+
+
+@boundscheck(False)
+@wraparound(False)
+cdef inline void _parallel_gain_thresholding(
+ raw_t[:, :, ::contiguous] in_raw,
+ int[:, ::contiguous] parallel_gain_indices,
+ data_t[:, :, :, ::contiguous] ccv_noise,
+ data_t sigma_high, data_t sigma_mid,
+ int frame, int cell, int ss, int fs,
+ data_t* data_ptr, gain_t* gain_ptr
+) noexcept nogil:
+ cdef int frame_high, frame_mid, frame_low
+ cdef data_t data_high, data_mid, data_low
+ cdef data_t threshold_high, threshold_mid
+
+ # Obtain indices to this pixel in each of three gain images.
+ frame_high = parallel_gain_indices[frame, 0]
+ frame_mid = parallel_gain_indices[frame, 1]
+ frame_low = parallel_gain_indices[frame, 2]
+
+ if (
+ ((in_raw[frame_high, ss, fs] & 0x3000) >> 12) == 0 and
+ ((in_raw[frame_mid, ss, fs] & 0x3000) >> 12) == 1 and
+ ((in_raw[frame_low, ss, fs] & 0x3000) >> 12) == 2
+ ):
+ # Verify that this pixel is recorded in the correct gain stage
+ # in each of the gain images. Memory cells in the transition
+ # regions between the gains can sometimes end up in the wrong
+ # one.
+
+ # Decode intensity in every gain stage.
+ data_high = <data_t>(in_raw[frame_high, ss, fs] & 0xFFF)
+ data_mid = <data_t>(in_raw[frame_mid, ss, fs] & 0xFFF)
+ data_low = <data_t>(in_raw[frame_low, ss, fs] & 0xFFF)
+
+ # Compute thresholds based on noise level.
+ threshold_high = 4096 - sigma_high * ccv_noise[cell, ss, fs, 0]
+ threshold_mid = 4096 - sigma_mid * ccv_noise[cell, ss, fs, 1]
+
+ # Pick the optimal gain stage for this pixel.
+ if data_mid > threshold_mid:
+ data_ptr[0] = data_low
+ gain_ptr[0] = 2
+ elif data_high > threshold_high:
+ data_ptr[0] = data_mid
+ gain_ptr[0] = 1
+ else:
+ data_ptr[0] = data_high
+ gain_ptr[0] = 0
+
+ else:
+ # Using an invalid gain stage triggers bad pixel masking later
+ # in the correction kernel.
+ gain_ptr[0] = 4