diff --git a/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb b/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb index 442f0d65e6b338125d96a4eb47fce961ca45b3f5..5ef0c5b28f63a598a4809547edaf703e47c9eeb5 100644 --- a/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb +++ b/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb @@ -89,6 +89,9 @@ "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", + "# Additional processing\n", + "count_lit_pixels = False # Count the number of pixels registering photons\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", @@ -259,6 +262,7 @@ " 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", + " corr_bools[\"count_lit_pixels\"] = count_lit_pixels\n", "\n", "# Many corrections don't apply to fixed gain mode; will explicitly disable later if detected\n", "disable_for_fixed_gain = [\n", @@ -532,7 +536,8 @@ "metadata": {}, "outputs": [], "source": [ - "if round_photons and photon_energy <= 0.0:\n", + "photon_energy_warn = None\n", + "if photon_energy <= 0.0:\n", " if use_xgm_device:\n", " # Try to obtain photon energy from XGM device.\n", " wavelength_data = dc[use_xgm_device, 'pulseEnergy.wavelengthUsed']\n", @@ -544,18 +549,57 @@ " photon_energy = (h * c / e) / (wavelength_data.as_single_value(rtol=1e-2) * 1e-6)\n", " print(f'Obtained photon energy {photon_energy:.3f} keV from {use_xgm_device}')\n", " except ValueError:\n", - " warning('XGM source available but photon energy varies greater than 1%, '\n", - " 'photon rounding disabled!')\n", - " round_photons = False\n", + " photon_energy_warn = 'XGM source available but photon energy varies greater than 1%'\n", " else:\n", - " warning('Neither explicit photon energy nor XGM device configured, photon rounding disabled!')\n", + " photon_energy_warn = 'Neither explicit photon energy nor XGM device configured'\n", + "\n", + "rounding_threshold_warn = None\n", + "if rounding_threshold <= .0 or 1. <= rounding_threshold:\n", + " rounding_threshold_warn = 'Round threshold is out of (0, 1) range. Use standard 0.5 value.'\n", + " rounding_threshold = 0.5\n", + "\n", + "if round_photons:\n", + " if photon_energy_warn:\n", + " warning(photon_energy_warn + ', photon rounding disabled!')\n", " round_photons = False\n", - "elif round_photons:\n", - " print(f'Photon energy for rounding: {photon_energy:.3f} keV')\n", + " else:\n", + " print(f'Photon energy for rounding: {photon_energy:.3f} keV')\n", + " if rounding_threshold_warn:\n", + " warning(rounding_threshold_warn)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "if round_photons:\n", + " data_units = 'photon'\n", + " litpx_threshold = 1.\n", + " if not xray_gain:\n", + " litpx_threshold *= 7.\n", + "else:\n", + " data_units = 'keV'\n", + " if photon_energy_warn:\n", + " warning(photon_energy_warn, '. Use 12 keV for lit-pixel counter threshold')\n", + " litpx_threshold = 12.\n", + " else:\n", + " litpx_threshold = photon_energy\n", + "\n", + " if rounding_threshold_warn:\n", + " warning(rounding_threshold_warn)\n", + "\n", + "if not xray_gain:\n", + " # convert photon energy to ADU (for AGIPD approx. 1 keV = 7 ADU)\n", + " # it looks that rounding to photons can be applied to data in ADU as well\n", + " litpx_threshold *= 7.\n", + " data_units = 'ADU'\n", "\n", - "if round_photons and (rounding_threshold <= .0 or 1. <= rounding_threshold):\n", - " warning('Round threshould is out of (0, 1) range. Use standard 0.5 value.')\n", - " rounding_threshold = 0.5" + "litpx_threshold *= rounding_threshold\n", + "\n", + "if count_lit_pixels:\n", + " print(f\"Count lit-pixels with signal above {litpx_threshold:.3g} {data_units}\")" ] }, { @@ -588,6 +632,7 @@ "agipd_corr.photon_energy = photon_energy\n", "agipd_corr.rounding_threshold = rounding_threshold\n", "agipd_corr.cs_mg_adjust = cs_mg_adjust\n", + "agipd_corr.litpx_threshold = litpx_threshold\n", "\n", "agipd_corr.compress_fields = compress_fields\n", "if recast_image_data:\n", @@ -906,6 +951,11 @@ " # Perform image-wise correction\"\n", " pool.starmap(agipd_corr.gain_correction, imagewise_chunks(img_counts))\n", " step_timer.done_step(\"Gain corrections\")\n", + " \n", + " # Peform additional processing\n", + " if count_lit_pixels:\n", + " pool.starmap(agipd_corr.litpixel_counter, imagewise_chunks(img_counts))\n", + " step_timer.done_step(\"Lit-pixel counting\")\n", "\n", " # Save corrected data\n", " pool.starmap(agipd_corr.write_file, [\n", diff --git a/src/cal_tools/agipdlib.py b/src/cal_tools/agipdlib.py index 05662d4fa126cca01f9f213668739e99d0b9c62c..7dd310093db6aacd4c4d5cc19fb74855d30cd640 100644 --- a/src/cal_tools/agipdlib.py +++ b/src/cal_tools/agipdlib.py @@ -582,6 +582,7 @@ class AgipdCorrections: self.photon_energy = 9.2 self.rounding_threshold = 0.5 self.cs_mg_adjust = 7e3 + self.litpx_threshold = 42 # Output parameters self.compress_fields = ['gain', 'mask'] @@ -614,7 +615,8 @@ 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', 'round_photons'] + 'low_medium_gap', 'round_photons', + 'count_lit_pixels'] if set(corr_bools).issubset(tot_corr_bools): self.corr_bools = corr_bools @@ -1495,10 +1497,49 @@ class AgipdCorrections: self.shared_dict[i]["valid_trains"] = sharedmem.empty(1024, dtype="u8") # noqa self.shared_dict[i]["nimg_in_trains"] = sharedmem.empty(1024, dtype="i8") # noqa + if self.corr_bools.get("count_lit_pixels"): + self.shared_dict[i]["litpx_counter"] = LitPixelCounter( + self.shared_dict[i], threshold=self.litpx_threshold) + if self.corr_bools.get("round_photons"): self.shared_hist_preround = sharedmem.empty(len(self.hist_bins_preround) - 1, dtype="i8") self.shared_hist_postround = sharedmem.empty(len(self.hist_bins_postround) - 1, dtype="i8") + def litpixel_counter(self, i_proc: int, first: int, last: int): + counter = self.shared_dict[i_proc].get("litpx_counter") + if counter: + counter.set_num_images(self.shared_dict[i_proc]["nImg"][0]) + counter.process(slice(first, last)) + + +class LitPixelCounter: + def __init__(self, data, threshold=0.8): + self.data = data + for name in ["data", "mask", "cellId", "pulseId", "trainId"]: + assert name in data + + self.image = data["data"] + self.mask = data["mask"] + + self.threshold = threshold + self.max_images = data["data"].shape[0] + + self.num_good_px = sharedmem.full(self.max_images, 0, int) + self.num_lit_px = sharedmem.full(self.max_images, 0, int) + + self.num_images = self.max_images + + def set_num_images(self, num_images): + self.num_images = num_images + + def process(self, chunk): + ix = range(*chunk.indices(self.num_images)) + for i in ix: + mask = self.mask[i] == 0 + self.num_lit_px[i] = np.sum( + self.image[i] > self.threshold, initial=0, where=mask) + self.num_good_px[i] = np.sum(mask) + def validate_selected_pulses( max_pulses: List[int], max_cells: int