From a521827a8d02afdc12374e36bdb58863fc294581 Mon Sep 17 00:00:00 2001
From: Philipp Schmidt <philipp.schmidt@xfel.eu>
Date: Mon, 11 Apr 2022 14:06:35 +0200
Subject: [PATCH] Add option to round intensity to absolute photon numbers in
AGIPD correct
---
.../AGIPD/AGIPD_Correct_and_Verify.ipynb | 41 +++++++++++++++++--
src/cal_tools/agipdlib.py | 31 +++++++++++++-
2 files changed, 68 insertions(+), 4 deletions(-)
diff --git a/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb b/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb
index 678a2b663..4684bdadf 100644
--- a/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb
+++ b/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb
@@ -40,9 +40,6 @@
"cal_db_timeout = 30000 # in milliseconds\n",
"creation_date_offset = \"00:00:00\" # add an offset to creation date, e.g. to get different constants\n",
"\n",
- "use_ppu_device = '' # Device ID for a pulse picker device to only process picked trains, empty string to disable\n",
- "ppu_train_offset = 0 # When using the pulse picker, offset between the PPU's sequence start and actually picked train\n",
- "\n",
"mem_cells = 0 # Number of memory cells used, set to 0 to automatically infer\n",
"bias_voltage = 0 # bias voltage, set to 0 to use stored value in slow data.\n",
"acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine\n",
@@ -84,10 +81,17 @@
"melt_snow = False # Identify (and optionally interpolate) 'snowy' pixels\n",
"mask_zero_std = False # Mask pixels with zero standard deviation across train\n",
"low_medium_gap = False # 5 sigma separation in thresholding between low and medium gain\n",
+ "round_photons = False # Round to absolute number of photons, only use with gain corrections\n",
+ "\n",
+ "# Optional auxiliary devices\n",
+ "use_ppu_device = '' # Device ID for a pulse picker device to only process picked trains, empty string to disable\n",
+ "ppu_train_offset = 0 # When using the pulse picker, offset between the PPU's sequence start and actually picked train\n",
"\n",
"use_litframe_device = '' # Device ID for a lit frame finder device to only process illuminated frames, empty string to disable\n",
"energy_threshold = -1000 # The low limit for the energy (uJ) exposed by frames subject to processing. If -1000, selection by pulse energy is disabled\n",
"\n",
+ "use_xgm_device = 'SA2_XTD1_XGM/XGM/DOOCS' # DoocsXGM device ID to obtain actual photon energy, operating condition else.\n",
+ "\n",
"# Output parameters\n",
"recast_image_data = '' # Cast data to a different dtype before saving\n",
"compress_fields = ['gain', 'mask'] # Datasets in image group to compress.\n",
@@ -218,6 +222,7 @@
" corr_bools[\"melt_snow\"] = melt_snow\n",
" corr_bools[\"mask_zero_std\"] = mask_zero_std\n",
" corr_bools[\"low_medium_gap\"] = low_medium_gap\n",
+ " corr_bools[\"round_photons\"] = round_photons\n",
"\n",
"# Many corrections don't apply to fixed gain mode; will explicitly disable later if detected\n",
"disable_for_fixed_gain = [\n",
@@ -462,6 +467,35 @@
"print(cell_sel.msg())"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "actual_photon_energy = None\n",
+ "\n",
+ "if use_xgm_device:\n",
+ " # Try to obtain photon energy from XGM device.\n",
+ " wavelength_data = dc[use_xgm_device, 'pulseEnergy.wavelengthUsed']\n",
+ " \n",
+ " try:\n",
+ " from scipy.constants import h, c, e\n",
+ " \n",
+ " # Read wavelength as a single value and convert to hv.\n",
+ " actual_photon_energy = (h * c / e) / (wavelength_data.as_single_value(rtol=1e-2) * 1e-6)\n",
+ " print(f'Obtained actual photon energy {actual_photon_energy:.3f} keV from {use_xgm_device}')\n",
+ " except ValueError:\n",
+ " if round_photons:\n",
+ " print('WARNING: XGM source available but actual photon energy varies greater than 1%, '\n",
+ " 'photon rounding disabled!')\n",
+ " round_photons = False\n",
+ "\n",
+ "if actual_photon_energy is None and round_photons:\n",
+ " print('WARNING: Using operating condition for actual photon energy in photon rounding mode, this is NOT reliable!')\n",
+ " actual_photon_energy = photon_energy"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
@@ -496,6 +530,7 @@
"agipd_corr.cm_n_itr = cm_n_itr\n",
"agipd_corr.noisy_adc_threshold = noisy_adc_threshold\n",
"agipd_corr.ff_gain = ff_gain\n",
+ "agipd_corr.actual_photon_energy = actual_photon_energy\n",
"\n",
"agipd_corr.compress_fields = compress_fields\n",
"if recast_image_data:\n",
diff --git a/src/cal_tools/agipdlib.py b/src/cal_tools/agipdlib.py
index 69f4d4aab..c86cd2d25 100644
--- a/src/cal_tools/agipdlib.py
+++ b/src/cal_tools/agipdlib.py
@@ -343,6 +343,7 @@ class AgipdCorrections:
self.hg_hard_threshold = 100
self.noisy_adc_threshold = 0.25
self.ff_gain = 1
+ self.actual_photon_energy = 9.2
# Output parameters
self.compress_fields = ['gain', 'mask']
@@ -367,7 +368,7 @@ class AgipdCorrections:
'zero_orange', 'force_hg_if_below',
'force_mg_if_below', 'mask_noisy_adc',
'melt_snow', 'common_mode', 'mask_zero_std',
- 'low_medium_gap']
+ 'low_medium_gap', 'round_photons']
if set(corr_bools).issubset(tot_corr_bools):
self.corr_bools = corr_bools
@@ -851,6 +852,34 @@ class AgipdCorrections:
msk[bidx] |= BadPixels.VALUE_OUT_OF_RANGE
del bidx
+ # Round keV-normalized intensity to photons.
+ if self.corr_bools.get("round_photons"):
+ data /= self.actual_photon_energy
+ np.round(data, out=data)
+
+ # This could also be done before and its mask inverted for
+ # rounding, but the performance difference is negligible.
+ bidx = data < 0
+ data[bidx] = 0
+ msk[bidx] |= BadPixels.VALUE_OUT_OF_RANGE
+ del bidx
+
+ if np.issubdtype(self.recast_image_fields.get('image'), np.integer):
+ # If the image data is meant to be recast to an integer
+ # type, make sure its values are within its bounds.
+
+ type_info = np.iinfo(self.recast_image_data['image'])
+
+ bidx = data < type_info.min
+ data[bidx] = 0
+ msk[bidx] |= BadPixels.VALUE_OUT_OF_RANGE
+ del bidx
+
+ bidx = data > type_info.max
+ data[bidx] = type_info.max
+ msk[bidx] |= BadPixels.VALUE_OUT_OF_RANGE
+ del bidx
+
# Mask entire ADC if they are noise above a threshold
# TODO: Needs clarification if needed,
# the returned arg is not used.
--
GitLab