From f9c883b5efecc04474b0883319ae02e705ebd4f6 Mon Sep 17 00:00:00 2001
From: Steffen Hauf <haufs@max-exfl059.desy.de>
Date: Tue, 13 Nov 2018 16:51:16 +0100
Subject: [PATCH] PEP8 compliance
---
cal_tools/cal_tools/agipdlib.py | 1385 +++++++++++++++++--------------
1 file changed, 784 insertions(+), 601 deletions(-)
diff --git a/cal_tools/cal_tools/agipdlib.py b/cal_tools/cal_tools/agipdlib.py
index 2aeef583d..bb17e4393 100644
--- a/cal_tools/cal_tools/agipdlib.py
+++ b/cal_tools/cal_tools/agipdlib.py
@@ -5,11 +5,11 @@ import h5py
import numpy as np
from scipy.signal import cwt, ricker
from sklearn.mixture import GaussianMixture
-from sklearn.preprocessing import StandardScaler
+from sklearn.preprocessing import StandardScaler
from cal_tools.enums import BadPixels
from cal_tools.tools import get_constant_from_db
-from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
+from iCalibrationDB import Constants, Conditions, Detectors
class SnowResolution(Enum):
@@ -17,7 +17,7 @@ class SnowResolution(Enum):
"""
NONE = "none"
INTERPOLATE = "interpolate"
-
+
class AgipdCorrections:
"""
@@ -29,12 +29,14 @@ class AgipdCorrections:
infile = h5py.File(filename, "r", driver="core")
outfile = h5py.File(filename_out, "w")
- agp_corr = AgipdCorrections(infile, outfile, max_cells, channel, max_pulses,
+ agp_corr = AgipdCorrections(infile, outfile, max_cells, channel,
+ max_pulses,
bins_gain_vs_signal, bins_signal_low_range,
bins_signal_high_range)
try:
- agp_corr.initialize(offset, rel_gain, rel_gain_offset, mask, noise, flatfield)
+ agp_corr.initialize(offset, rel_gain, rel_gain_offset, mask,
+ noise, flatfield)
except IOError:
return
@@ -43,16 +45,18 @@ class AgipdCorrections:
hists, edges = agp_corr.get_histograms()
- """
+ """
+
def __init__(self, infile, outfile, max_cells, channel, max_pulses,
- bins_gain_vs_signal, bins_signal_low_range, bins_signal_high_range,
+ bins_gain_vs_signal, bins_signal_low_range,
+ bins_signal_high_range,
bins_dig_gain_vs_signal, raw_fmt_version=2, chunk_size=512,
h5_data_path="INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/",
h5_index_path="INDEX/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/",
do_rel_gain=True, chunk_size_idim=512):
"""
Initialize an AgipdCorrections Class
-
+
:param infile: to be corrected h5py input file
:param outfile: writeable h5py output file
:param max_cell: maximum number of memory cells to handle, e.g. if
@@ -60,14 +64,19 @@ class AgipdCorrections:
:param channel: module/channel to correct
:param max_pulses: maximum pulse id to consider for preview histograms
:param bins_gain_vs_signal: number of bins for gain vs signal histogram
- :param bins_signal_low_range: number of bins for the low signal range histogram
- :param bins_signal_high_range: number of bins for the high signal range histogram
+ :param bins_signal_low_range: number of bins for the low signal
+ range histogram
+ :param bins_signal_high_range: number of bins for the high signal
+ range histogram
:param raw_fmt_version: raw file format version to use
:param chunk_size: images per chunk to compute upon
- :param h5_data_path: path in HDF5 file which is prefixed to the image/data section
- :param h5_index_path: path in HDF5 file which is prefixed to the index section
+ :param h5_data_path: path in HDF5 file which is prefixed to the
+ image/data section
+ :param h5_index_path: path in HDF5 file which is prefixed to the
+ index section
:param do_rel_gain: do relative gain corrections
- :param chunk_size_idim: chunking size on image dimension when writing data out
+ :param chunk_size_idim: chunking size on image dimension when
+ writing data out
"""
self.agipd_base = h5_data_path.format(channel)
self.idx_base = h5_index_path.format(channel)
@@ -98,161 +107,193 @@ class AgipdCorrections:
self.melt_snow = SnowResolution.NONE
self.match_asics = True
self.chunk_size_idim = chunk_size_idim
-
+ self.dohigh = 0
+
def get_iteration_range(self):
- """Returns a range expression over which to iterate in chunks
+ """Returns a range expression over which to iterate in chunks
"""
- return np.array_split(self.firange, self.firange.size//self.chunksize)
-
- def initialize(self, offset, rel_gain, mask, noise, thresholds, base_offset, swap_axis=True):
+ return np.array_split(self.firange,
+ self.firange.size // self.chunksize)
+
+ def initialize(self, offset, rel_gain, mask, noise, thresholds,
+ base_offset, swap_axis=True):
""" Initialize the calibration constants and the output data file.
-
- Any data that is not touched by the corrections is copied into the output
- file at this point. Also data validity tests are performed. This functions
- must be called before `correct_lpd` is executed.
-
+
+ Any data that is not touched by the corrections is copied into the
+ output file at this point. Also data validity tests are performed.
+ This functions must be called before `correct_lpd` is executed.
+
:param offset: offset map to use for corrections
:param rel_gain: relative gain map to use for corrections
:param mask: bad pixel mask to use for corrections
:param noise: noise map to use for corrections
:param thresholds: thresholds for gain determination
- :param base_offset: base offsets from PC and FF data: should be a tuple:
- PC offset for high gain, PC offset for medium gain, FF offset for high gain
- :param swap_axis: set to True if using data from the calibration database.
-
- Note that this function may be called twice, when initializing partially from DB and file.
+ :param base_offset: base offsets from PC and FF data: should be a
+ tuple:
+ PC offset for high gain, PC offset for medium gain, FF offset
+ for high gain
+ :param swap_axis: set to True if using data from the calibration
+ database.
+
+ Note that this function may be called twice, when initializing
+ partially from DB and file.
Constants not to be set by the initial call should be set to `None`.
"""
-
- # swap the axes such that memory cell dimension is first - the usual way of operation,
- # and the one expected by the rest of this class.
+
+ # swap the axes such that memory cell dimension is first - the usual
+ # way of operation, and the one expected by the rest of this class.
if swap_axis:
if offset is not None:
- offset = np.ascontiguousarray(np.moveaxis(np.moveaxis(offset, 2, 0), 2, 1))
+ mvd = np.moveaxis(np.moveaxis(offset, 2, 0), 2, 1)
+ offset = np.ascontiguousarray(mvd)
if rel_gain is not None:
- rel_gain = np.ascontiguousarray(np.moveaxis(np.moveaxis(rel_gain, 2, 0), 2, 1))
+ mvd = np.moveaxis(np.moveaxis(rel_gain, 2, 0), 2, 1)
+ rel_gain = np.ascontiguousarray(mvd)
if mask is not None:
- mask = np.ascontiguousarray(np.moveaxis(np.moveaxis(mask, 2, 0), 2, 1))
+ mvd = np.moveaxis(np.moveaxis(mask, 2, 0), 2, 1)
+ mask = np.ascontiguousarray(mvd)
if noise is not None:
- noise = np.ascontiguousarray(np.moveaxis(np.moveaxis(noise, 2, 0), 2, 1))
+ mvd = np.moveaxis(np.moveaxis(noise, 2, 0), 2, 1)
+ noise = np.ascontiguousarray(mvd)
if base_offset is not None:
for i in range(len(base_offset)):
if base_offset[i] is not None:
- base_offset[i] = np.ascontiguousarray(np.moveaxis(np.moveaxis(base_offset[i], 2, 0), 2, 1))
+ mvd = np.moveaxis(base_offset[i], 2, 0)
+ mvd = np.moveaxis(mvd, 2, 1)
+ base_offset[i] = np.ascontiguousarray(mvd)
if thresholds is not None:
- thresholds = np.ascontiguousarray(np.moveaxis(np.moveaxis(thresholds, 2, 0), 2, 1))
-
+ mvd = np.moveaxis(np.moveaxis(thresholds, 2, 0), 2, 1)
+ thresholds = np.ascontiguousarray(mvd)
+
if base_offset is not None:
self.base_offset = base_offset
if offset is not None:
self.offset = offset
-
+
if rel_gain is not None:
- self.rel_gain = 1./rel_gain
-
- # if a mask already exists we OR the mask to maintain what we already have
+ self.rel_gain = 1. / rel_gain
+
+ # if a mask already exists we OR the mask to maintain what we
+ # already have
if mask is not None:
if not hasattr(self, "mask"):
self.mask = mask
else:
- self.mask |= mask[:self.mask.shape[0],...]
-
+ self.mask |= mask[:self.mask.shape[0], ...]
+
if noise is not None:
self.noise = noise
-
+
if thresholds is not None:
self.thresholds = thresholds
-
+
if self.base_offset is not None and self.offset is not None:
# calculate a delta offset from base offsets and normal offsets
dohigh = np.zeros(self.offset.shape[:-1], np.float32)
for i in range(8):
for j in range(2):
- asic_m = np.nanmedian(self.offset[:, i*64:(i+1)*64, j*64:(j+1)*64, 0], axis=(1,2))
- global_m = np.nanmedian(self.offset[:, ..., 0], axis=(1,2))
- dohigh[:, i*64:(i+1)*64, j*64:(j+1)*64] = (asic_m-global_m)[:,None,None]
+ asico = self.offset[:, i * 64:(i + 1) * 64,
+ j * 64:(j + 1) * 64, 0]
+ asic_m = np.nanmedian(asico, axis=(1, 2))
+ global_m = np.nanmedian(self.offset[:, ..., 0],
+ axis=(1, 2))
+ doa = (asic_m - global_m)[:, None, None]
+ dohigh[:, i * 64:(i + 1) * 64, j * 64:(j + 1) * 64] = doa
+
self.dohigh = dohigh
-
- dopc = self.base_offset[1]-self.base_offset[0]
- dooff = self.offset[...,1] - self.offset[...,0]
+
+ dopc = self.base_offset[1] - self.base_offset[0]
+ dooff = self.offset[..., 1] - self.offset[..., 0]
doff = dopc - dooff
for i in range(8):
for j in range(2):
- med = np.nanmedian(doff[:, i*64:(i+1)*64, j*64:(j+1)*64], axis=(1,2))
- doff[:, i*64:(i+1)*64, j*64:(j+1)*64] = med[:, None, None]
- self.offset[...,1] += doff
-
+ ado = doff[:, i * 64:(i + 1) * 64, j * 64:(j + 1) * 64]
+ med = np.nanmedian(ado, axis=(1, 2))[:, None, None]
+ doff[:, i * 64:(i + 1) * 64, j * 64:(j + 1) * 64] = med
+
+ self.offset[..., 1] += doff
+
# this is the first call, usually from the database
# we also generate file structures here, as the second call is optional
if not self.initialized:
- self.median_noise = np.nanmedian(self.noise[0,...])
- self.max_cells = np.min([self.offset.shape[0], self.mask.shape[0], self.max_cells])
- self.gen_valid_range()
+ self.median_noise = np.nanmedian(self.noise[0, ...])
+ cops = [self.offset.shape[0], self.mask.shape[0], self.max_cells]
+ self.max_cells = np.min(cops)
+ self.gen_valid_range()
self.copy_and_sanitize_non_cal_data()
self.create_output_datasets()
self.initialized = True
- print("Offet medians are {}".format(np.nanmedian(self.offset, axis=(0,1,2))))
+ print("Offet medians are {}".format(
+ np.nanmedian(self.offset, axis=(0, 1, 2))))
if hasattr(self, "rel_gain"):
- print("Gain medians are {}".format(np.nanmedian(self.rel_gain, axis=(0,1,2))))
- print("Threshold medians are {}".format(np.nanmedian(self.thresholds, axis=(0,1,2))))
+ print("Gain medians are {}".format(
+ np.nanmedian(self.rel_gain, axis=(0, 1, 2))))
+ print("Threshold medians are {}".format(
+ np.nanmedian(self.thresholds, axis=(0, 1, 2))))
else:
print("After reading from file: ")
- print("Offet medians are {}".format(np.nanmedian(self.offset, axis=(0,1,2))))
- print("Gain medians are {}".format(np.nanmedian(self.rel_gain, axis=(0,1,2))))
- print("Threshold medians are {}".format(np.nanmedian(self.thresholds, axis=(0,1,2))))
-
+ print("Offet medians are {}".format(
+ np.nanmedian(self.offset, axis=(0, 1, 2))))
+ print("Gain medians are {}".format(
+ np.nanmedian(self.rel_gain, axis=(0, 1, 2))))
+ print("Threshold medians are {}".format(
+ np.nanmedian(self.thresholds, axis=(0, 1, 2))))
+
def split_gain(self, d):
""" Split gain information off AGIPD Data
"""
- return d[:,0,...], d[:,1,...]
-
+ return d[:, 0, ...], d[:, 1, ...]
+
def baseline_correct_via_noise(self, d, noise, g):
- """ Correct baseline shifts by evaluating the position of the noise peak
-
+ """ Correct baseline shifts by evaluating position of the noise peak
+
:param: d: the data to correct, should be a single image
:param noise: the mean noise for the cell id of `d`
:param g: gain array matching `d` array
-
- Correction is done by identifying the left-most (significant) peak in the
- histogram of `d` and shifting it to be centered on 0. This is done via
- a continous wavelet transform (CWT), using a Ricker wavelet.
-
- Only high gain data is evaulated, and data larger than 50 ADU, equivalent of
- rougly a 9 keV photon is ignored.
-
- Two passes are executed: the first shift is accurate to 10 ADU and will catch
- baseline shifts smaller then -2000 ADU. CWT is evaluated for peaks of widths
- one, three and five time the noise. The baseline is then shifted
- by the position of the summed CWT maximum.
-
+
+ Correction is done by identifying the left-most (significant) peak
+ in the histogram of `d` and shifting it to be centered on 0.
+ This is done via a continous wavelet transform (CWT), using a Ricker
+ wavelet.
+
+ Only high gain data is evaulated, and data larger than 50 ADU,
+ equivalent of roughly a 9 keV photon is ignored.
+
+ Two passes are executed: the first shift is accurate to 10 ADU and
+ will catch baseline shifts smaller then -2000 ADU. CWT is evaluated
+ for peaks of widths one, three and five time the noise.
+ The baseline is then shifted by the position of the summed CWT maximum.
+
In a second pass histogram is evaluated within a range of
+- 5*noise of the initial shift value, for peak of width noise.
-
+
Baseline shifts larger than the maximum threshold or positive shifts
larger 10 are discarded, and the original data is return, otherwise
the shift corrected data in returned.
-
+
"""
- # we work on a copy to be able to filter values by setting them to np.nan
+ # we work on a copy to be able to filter values by setting them to
+ # np.nan
dd = copy.copy(d)
dd[g != 0] = np.nan # only high gain data
dd[dd > 50] = np.nan # only noise data
h, e = np.histogram(dd, bins=210, range=(-2000, 100))
- c = (e[1:]+e[:-1])/2
+ c = (e[1:] + e[:-1]) / 2
try:
- cwtmatr = cwt(h, ricker, [noise, 3.*noise, 5.*noise])
+ cwtmatr = cwt(h, ricker, [noise, 3. * noise, 5. * noise])
except:
return d
cwtall = np.sum(cwtmatr, axis=0)
marg = np.argmax(cwtall)
pc = c[marg]
- high_res_range = (dd > pc - 5*noise) & (dd < pc + 5*noise)
+ high_res_range = (dd > pc - 5 * noise) & (dd < pc + 5 * noise)
dd[~high_res_range] = np.nan
- h, e = np.histogram(dd, bins=100, range=(pc - 5*noise, pc + 5*noise))
- c = (e[1:]+e[:-1])/2
+ h, e = np.histogram(dd, bins=100,
+ range=(pc - 5 * noise, pc + 5 * noise))
+ c = (e[1:] + e[:-1]) / 2
try:
cwtmatr = cwt(h, ricker, [noise, ])
except:
@@ -262,89 +303,92 @@ class AgipdCorrections:
# too large shift to be easily decernable via noise
if pc > 10 or pc < -self.baseline_corr_noise_threshold:
return d
- return d-pc
+ return d - pc
def correct_baseline_via_hist(self, d, pcm, g):
- """ Correct baseline shifts by matching edges of high and medium gain histogram
-
+ """ Correct baseline shifts by matching edges of high and medium
+ gain histogram
+
:param d: single image to correct
:param pcm: relative gain slope for medium gain of each pixel in `d`
:param g: gain values for `d`
-
- As a preparation the median value of medium gain pixels is shifted in
+
+ As a preparation the median value of medium gain pixels is shifted in
increments of 50 ADU as long as it is negative.
-
+
Correction is then performed by evaluating histograms for high gain
- and medium gain pixels in `d`. The right-most significant bin of the high
- gain histogram is then shifted such that it coincides with the left-most
- significant bin of the medium gain histogram. Significant here means that
- bin counts are at least 10% of the average bin count of histogram for all
- bins larger 10 counts. This initial evaluation uses histograms in range
- between -10000 and +10000 ADU with a resultion of 100 ADU per bin. The
- shift required to match both histogram borders is an initial estimate for the
- baseline shift.
-
+ and medium gain pixels in `d`. The right-most significant bin of the
+ high gain histogram is then shifted such that it coincides with the
+ left-most significant bin of the medium gain histogram. Significant
+ here means that bin counts are at least 10% of the average bin count
+ of histogram for all bins larger 10 counts. This initial evaluation
+ uses histograms in range between -10000 and +10000 ADU with a
+ resultion of 100 ADU per bin. The shift required to match both
+ histogram borders is an initial estimate for the baseline shift.
+
Finally, the mean bin count of the five outermost bins of the high and
- medium gain histograms is compared. The baseline is shifted by increments
- of 1 ADU until they are within 20% of each other.
-
- From this point on additional shifts are performed while the deviation stays
- within 30% of each other. The final shift is then evaluated as the point of minimal
- deviation of these values.
-
- A maximum iteration count of 200 is imposed on the procedure. If no convergence
- is reached within this limit, the original data is returned, otherwise the
- baseline corrected data is returned.
+ medium gain histograms is compared. The baseline is shifted by
+ increments of 1 ADU until they are within 20% of each other.
+
+ From this point on additional shifts are performed while the
+ deviation stays within 30% of each other. The final shift is then
+ evaluated as the point of minimal deviation of these values.
+
+ A maximum iteration count of 200 is imposed on the procedure. If no
+ convergence is reached within this limit, the original data is
+ returned, otherwise the baseline corrected data is returned.
"""
dd = copy.copy(d)
# shift while the medium gain produces negative values on average
pc = 0
it = 0
max_it = 100
- while np.nanmedian(dd[g==1]-pc) < 0:
+ while np.nanmedian(dd[g == 1] - pc) < 0:
pc -= 50
if it > max_it:
return d
it += 1
- def min_hist_distance(pc, bins=100, ran=(-10000, 10000), dec=20, minbin=10):
- opc = pc
- hh, e = np.histogram(dd[g==0]-pc, bins=bins, range=ran)
- hm, e = np.histogram((dd[g==1]-pc)*pcm[g==1], bins=bins, range=ran)
+ def min_hist_distance(pc, bins=100, ran=(-10000, 10000), dec=20,
+ minbin=10):
+ hh, e = np.histogram(dd[g == 0] - pc, bins=bins, range=ran)
+ hm, e = np.histogram((dd[g == 1] - pc) * pcm[g == 1], bins=bins,
+ range=ran)
# right most significant value of high gain histogram
hhm = np.nanmean(hh[hh > minbin])
- rmh = hh.size-np.argmax((hh > 0.1*hhm)[::-1])
+ rmh = hh.size - np.argmax((hh > 0.1 * hhm)[::-1])
# left most significant value of high gain histogram
hmm = np.nanmean(hm[hm > minbin])
- lmm = np.argmax((hm > 0.1*hmm))
- med_pcm = np.nanmedian(pcm[g==1])
- eh = e[rmh]
+ lmm = np.argmax((hm > 0.1 * hmm))
+ med_pcm = np.nanmedian(pcm[g == 1])
+ eh = e[rmh]
el = e[lmm]
- pc += (el - eh)/(med_pcm -1)
+ pc += (el - eh) / (med_pcm - 1)
return pc
# set a pc one halfstep higher than initial guesstimate
pc += 50
pc = min_hist_distance(pc)
-
+
# now start rolling back until we have similar signal levels
def comp_signal_level(pc, minbin=1):
- hh, e = np.histogram(dd[g==0]-pc, bins=100, range=(5000, 7000))
- hm, e = np.histogram((dd[g==1]-pc)*pcm[g==1], bins=100, range=(5000, 7000))
+ hh, e = np.histogram(dd[g == 0] - pc, bins=100, range=(5000, 7000))
+ hm, e = np.histogram((dd[g == 1] - pc) * pcm[g == 1], bins=100,
+ range=(5000, 7000))
# right most significant value of high gain histogram
hhm = np.nanmean(hh[hh > minbin])
- rmh = hh.size-np.argmax((hh > 0.1*hhm)[::-1])
+ rmh = hh.size - np.argmax((hh > 0.1 * hhm)[::-1])
# left most significant value of high gain histogram
hmm = np.nanmean(hm[hm > minbin])
- lmm = np.argmax((hm > 0.1*hmm))
+ lmm = np.argmax((hm > 0.1 * hmm))
reg = 5
- ret = np.abs(np.sum(hh[rmh-reg:rmh])-np.sum(hm[lmm:lmm+reg]))
- ret /= np.sum(hh[rmh-reg:rmh])
+ ret = np.abs(np.sum(hh[rmh - reg:rmh]) - np.sum(hm[lmm:lmm + reg]))
+ ret /= np.sum(hh[rmh - reg:rmh])
return ret
it = 0
@@ -365,7 +409,7 @@ class AgipdCorrections:
slevs = []
slev = comp_signal_level(pc)
slevs.append(slev)
- while slev < 0.3:
+ while slev < 0.3:
pc += 1
slev = comp_signal_level(pc)
slevs.append(slev)
@@ -374,61 +418,64 @@ class AgipdCorrections:
pc = pcs[np.argmin(slevs)]
- return d-pc
-
+ return d - pc
+
def correct_baseline_via_hist_asic(self, d, g):
""" Correct diagonal falloff on ASICs by matching corner histograms
-
+
:param d: single image data
:param g: gain values for `d`
-
- Corrections are performed for the top and bottom row of ASICs seperately.
-
- In easch row, starting from the beam-hole closes ASIC, the histogram of a
- square region of size 8x8 pixels is evaluted for its maximum bin count
- location (only medium gain values), and compared to a histogram produced
- from a neighbouring region on the next ASIC. Double sized pixels are avoided.
-
- The reasoning is that the histograms should not dramatically differ for
- continously distributed photon events. Each ASIC is then shifted such that
- histograms much and the corrected data is returned.
-
+
+ Corrections are performed for the top and bottom row of ASICs
+ seperately.
+
+ In easch row, starting from the beam-hole closes ASIC, the histogram
+ of a square region of size 8x8 pixels is evaluted for its maximum
+ bin count location (only medium gain values), and compared to a
+ histogram produced from a neighbouring region on the next ASIC.
+ Double sized pixels are avoided.
+
+ The reasoning is that the histograms should not dramatically differ
+ for continuously distributed photon events. Each ASIC is then shifted
+ such that histograms much and the corrected data is returned.
+
Warning: this feature is very experimental!
"""
-
+
rsize = 8
- def hist_ends(dm, bins=5000, ran=(-5000, 5000), minbin=4/(16/rsize)):
- hm, e = np.histogram(dm, bins=bins, range=ran)
+ def hist_ends(dm, bins=5000, ran=(-5000, 5000),
+ minbin=4 / (16 / rsize)):
+ hm, e = np.histogram(dm, bins=bins, range=ran)
# left most significant value of medium gain histogram
- #hmm = np.nanmean(hm[hm > minbin])
+ # hmm = np.nanmean(hm[hm > minbin])
lmm = np.argmax((hm > minbin))
- return e[lmm], np.sum(hm)
-
+ return e[lmm], np.sum(hm)
for i in range(1, 8):
# upper asics
- casic = np.concatenate((d[i*64-rsize:i*64-1,64:64+rsize].flatten(),
- d[i*64+1:i*64+rsize,64-rsize:64].flatten()))
-
- casic_g = np.concatenate((g[i*64-rsize:i*64-1,64:64+rsize].flatten(),
- g[i*64+1:i*64+rsize,64-rsize:64].flatten()))
+ casic = np.concatenate(
+ (d[i * 64 - rsize:i * 64 - 1, 64:64 + rsize].flatten(),
+ d[i * 64 + 1:i * 64 + rsize, 64 - rsize:64].flatten()))
+ casic_g = np.concatenate(
+ (g[i * 64 - rsize:i * 64 - 1, 64:64 + rsize].flatten(),
+ g[i * 64 + 1:i * 64 + rsize, 64 - rsize:64].flatten()))
idxm = casic_g == 1
cme, cms = hist_ends(casic[idxm])
- asic = d[i*64+1:i*64+rsize,64:64+rsize].flatten()
- asic_g = g[i*64+1:i*64+rsize,64:64+rsize].flatten()
+ asic = d[i * 64 + 1:i * 64 + rsize, 64:64 + rsize].flatten()
+ asic_g = g[i * 64 + 1:i * 64 + rsize, 64:64 + rsize].flatten()
idxm = asic_g == 1
- me, ms = hist_ends( asic[idxm])
+ me, ms = hist_ends(asic[idxm])
- if cms > rsize*rsize/10 and ms > rsize*rsize/10:
+ if cms > rsize * rsize / 10 and ms > rsize * rsize / 10:
pc = cme - me
else:
pc = 0
@@ -436,28 +483,32 @@ class AgipdCorrections:
if np.abs(pc) > 500:
# somthing when wrong
pc = 0
- d[i*64:(i+1)*64,64:] += pc
-
+ d[i * 64:(i + 1) * 64, 64:] += pc
for i in range(0, 7):
# lower asics
- casic = np.concatenate((d[(i+1)*64-rsize:(i+1)*64-1,64:64+rsize].flatten(),
- d[(i+1)*64+1:(i+1)*64+rsize,64-rsize:64].flatten()))
-
- casic_g = np.concatenate((g[(i+1)*64-rsize:(i+1)*64-1,64:64+rsize].flatten(),
- g[(i+1)*64+1:(i+1)*64+rsize,64-rsize:64].flatten()))
+ casic = np.concatenate((d[(i + 1) * 64 - rsize:(i + 1) * 64 - 1,
+ 64:64 + rsize].flatten(),
+ d[(i + 1) * 64 + 1:(i + 1) * 64 + rsize,
+ 64 - rsize:64].flatten()))
+ casic_g = np.concatenate((g[(i + 1) * 64 - rsize:(i + 1) * 64 - 1,
+ 64:64 + rsize].flatten(),
+ g[(i + 1) * 64 + 1:(i + 1) * 64 + rsize,
+ 64 - rsize:64].flatten()))
idxm = casic_g == 1
cme, cms = hist_ends(casic[idxm])
- asic = d[(i+1)*64-rsize:(i+1)*64-1,64-rsize:64].flatten()
- asic_g = g[(i+1)*64-rsize:(i+1)*64-1,64-rsize:64].flatten()
+ asic = d[(i + 1) * 64 - rsize:(i + 1) * 64 - 1,
+ 64 - rsize:64].flatten()
+ asic_g = g[(i + 1) * 64 - rsize:(i + 1) * 64 - 1,
+ 64 - rsize:64].flatten()
idxm = asic_g == 1
- me, ms = hist_ends( asic[idxm])
+ me, ms = hist_ends(asic[idxm])
- if cms > rsize*rsize/10 and ms > rsize*rsize/10:
+ if cms > rsize * rsize / 10 and ms > rsize * rsize / 10:
pc = cme - me
else:
pc = 0
@@ -465,84 +516,97 @@ class AgipdCorrections:
if np.abs(pc) > 500:
# somthing went wrong
pc = 0
- d[i*64:(i+1)*64,:64] += pc
+ d[i * 64:(i + 1) * 64, :64] += pc
return d
-
+
def match_asic_borders(self, d, chunk=8):
""" Match border pixels of the two ASIC rows to the same median value
-
+
:param d: single image data
- :param chunk: number of pixels on each ASIC boundary to generate mean values for
-
- Each ASIC has 64//chunk mean value calculated along its two border pixel rows. The
- deviaiton of chunk pairs is then evaluated and the upper/lower ASICs are shifted
- by the mean deviation for the right/left hemisphere of the detector.
-
+ :param chunk: number of pixels on each ASIC boundary to generate
+ mean values for
+
+ Each ASIC has 64//chunk mean value calculated along its two border
+ pixel rows. The deviation of chunk pairs is then evaluated and the
+ upper/lower ASICs are shifted by the mean deviation for the
+ right/left hemisphere of the detector.
+
The corrected data is returned.
- """
+ """
for i in range(8):
- these_asics = d[:, i*64:(i+1)*64,:]
+ these_asics = d[:, i * 64:(i + 1) * 64, :]
diffs = np.zeros((d.shape[0], 8))
- for k in range(64//chunk):
-
- lowerasic = np.nanmedian(these_asics[:,k*chunk:(k+1)*chunk,62:64], axis=(1,2))
- upperasic = np.nanmedian(these_asics[:,k*chunk:(k+1)*chunk,64:66], axis=(1,2))
- diffs[:, k] = lowerasic-upperasic if self.channel < 8 else upperasic-lowerasic
+ for k in range(64 // chunk):
+ ldat = these_asics[:, k * chunk:(k + 1) * chunk, 62:64]
+ lowerasic = np.nanmedian(ldat, axis=(1, 2))
+ udat = these_asics[:, k * chunk:(k + 1) * chunk, 64:66]
+ upperasic = np.nanmedian(udat, axis=(1, 2))
+ low_up = lowerasic - upperasic
+ up_low = upperasic - lowerasic
+ diff = low_up if self.channel < 8 else up_low
+ diffs[:, k] = diff
+
+ mn = np.nanmean(diffs, axis=1)[:, None, None] / 2
if self.channel < 8:
- d[:, i*64:(i+1)*64,64:] += np.nanmean(diffs, axis=1)[:,None,None]/2
- d[:, i*64:(i+1)*64,:64] -= np.nanmean(diffs, axis=1)[:,None,None]/2
+ d[:, i * 64:(i + 1) * 64, 64:] += mn
+ d[:, i * 64:(i + 1) * 64, :64] -= mn
else:
- d[:, i*64:(i+1)*64,:64] += np.nanmean(diffs, axis=1)[:,None,None]/2
- d[:, i*64:(i+1)*64,64:] -= np.nanmean(diffs, axis=1)[:,None,None]/2
+ d[:, i * 64:(i + 1) * 64, :64] += mn
+ d[:, i * 64:(i + 1) * 64, 64:] -= mn
return d
-
+
def melt_snowy_pixels(self, raw, im, gain, rgain, resolution=None):
""" Identify (and optionally interpolate) 'snowy' pixels
-
- :param raw: raw data
+
+ :param raw: raw data
:param im: offet and relativegain corrected data:
:param gain: gain values for `raw`
- :param rgain: raw gain data, scaled by the threshold for high-to-medium gain
- switching
- :param resolution: resolution strategy, should of enum-type `SnowResolution`
-
- Snowy pixels are pixels which are already identified as pixels in the medium
- gain stage by their gain values, but which have transitional image values between
- the largest high gain value and the smalles medium gain value. As these value
- may be encountered again in the context of medium gain, they are ambigigous and
- cannot readily be identified by simple thresholding alone.
-
- It is attempted to identify these snowy pixels by using a Gaussian mixture
- clustering algorithm acting on multispectral imput. Positions in the cluster are
- identified as follows:
-
+ :param rgain: raw gain data, scaled by the threshold for
+ high-to-medium gain switching
+ :param resolution: resolution strategy, should of enum-type
+ `SnowResolution`
+
+ Snowy pixels are pixels which are already identified as pixels in
+ the medium gain stage by their gain values, but which have
+ transitional image values between the largest high gain value and
+ the smalles medium gain value. As these value may be encountered again
+ in the context of medium gain, they are ambigigous and cannot
+ readily be identified by simple thresholding alone.
+
+ It is attempted to identify these snowy pixels by using a Gaussian
+ mixture clustering algorithm acting on multispectral imput.
+ Positions in the cluster are identified as follows:
+
x: abs(p_r-p_bg)*p_rg**2
y: p_r
-
- where p_r is a given pixel raw value, p_bg is the mean value of the 8 direct
- neighbor pixels, and p_rg is the raw gain value of the pixel.
-
- Note that these cluster params are purely emprically derived, and do not have
- connection to ASIC design etc.
-
- Only pixels in medium gain are evaluated, and evaluation is image-wise, but
- subdivided further into ASICs.
-
- The so-created graph [[x_0, x_1, ..., x_n], [y_0, y_1, ..., y_n]] is scaled
- using a `sklearn.StandardScaler` and input into a two-component
- `GaussianMixture` clustering algorithm. This results in a set of labels,
- identifying pixels as likely snowy, or not. However, for a given image we
- do not know which label is which from the output. Hence, labels are
- differentiated under the assumption that snowy pixel occur be out-of-context,
- i.e. their surround pixels more likely at lower signal values, than if the
- pixel is in a region were gain switching led to a large value.
-
- Depending on the resolution strategy so-identified pixels are either set
- to `np.nan` or to the interpolated value of the direct (high-gain) neighbors.
-
- The corrected data is returned alonside an updated gain mask and bad pixel
+
+ where p_r is a given pixel raw value, p_bg is the mean value of the
+ 8 direct neighbor pixels, and p_rg is the raw gain value of the pixel.
+
+ Note that these cluster params are purely emprically derived,
+ and do not have connection to ASIC design etc.
+
+ Only pixels in medium gain are evaluated, and evaluation is
+ image-wise, but subdivided further into ASICs.
+
+ The so-created graph [[x_0, x_1, ..., x_n], [y_0, y_1, ..., y_n]] is
+ scaled using a `sklearn.StandardScaler` and input into a two-component
+ `GaussianMixture` clustering algorithm. This results in a set of
+ labels, identifying pixels as likely snowy, or not. However,
+ for a given image we do not know which label is which from the
+ output. Hence, labels are differentiated under the assumption that
+ snowy pixel occur be out-of-context, i.e. their surround pixels more
+ likely at lower signal values, than if the pixel is in a region were
+ gain switching led to a large value.
+
+ Depending on the resolution strategy so-identified pixels are either
+ set to `np.nan` or to the interpolated value of the direct (high-gain)
+ neighbors.
+
+ The corrected data is returned alonside an updated gain mask and bad
+ pixel
mask.
"""
snow_mask = np.zeros(im.shape, np.uint32)
@@ -550,22 +614,26 @@ class AgipdCorrections:
for i in range(8):
for j in range(2):
- fidx = gain[k, i*64:(i+1)*64,j*64:(j+1)*64] == 1
- fidx_h = gain[k, i*64:(i+1)*64,j*64:(j+1)*64] == 0
-
- # need at least two pixels in medium gain to be able to cluster in two groups
+ fidx = gain[k, i * 64:(i + 1) * 64,
+ j * 64:(j + 1) * 64] == 1
+ fidx_h = gain[k, i * 64:(i + 1) * 64,
+ j * 64:(j + 1) * 64] == 0
+
+ # need at least two pixels in medium gain to be able to
+ # cluster in two groups
if np.count_nonzero(fidx) < 2:
continue
# data for a given asic
- asic = raw[k, i*64:(i+1)*64,j*64:(j+1)*64]
- asic_g = gain[k, i*64:(i+1)*64,j*64:(j+1)*64]
- asic_r = rgain[k, i*64:(i+1)*64,j*64:(j+1)*64]
-
+ asic = raw[k, i * 64:(i + 1) * 64, j * 64:(j + 1) * 64]
+ asic_g = gain[k, i * 64:(i + 1) * 64, j * 64:(j + 1) * 64]
+ asic_r = rgain[k, i * 64:(i + 1) * 64, j * 64:(j + 1) * 64]
+
asic_h = copy.copy(asic)
asic_h[~fidx_h] = np.nan
-
- # calculate the background of by averaging the 8 direct neighbours of each pixel
+
+ # calculate the background of by averaging the 8 direct
+ # neighbours of each pixel
rl = np.roll(asic, 1, axis=0)
rr = np.roll(asic, -1, axis=0)
ru = np.roll(asic, 1, axis=1)
@@ -574,9 +642,10 @@ class AgipdCorrections:
rld = np.roll(rl, -1, axis=1)
rru = np.roll(rr, 1, axis=1)
rrd = np.roll(rr, -1, axis=1)
- bg = (rl + rr + ru + rd + rlu + rld + rru + rrd)/8
+ bg = (rl + rr + ru + rd + rlu + rld + rru + rrd) / 8
- # calculate the background of by averaging the 8 direct neighbours of each pixel
+ # calculate the background of by averaging the 8 direct
+ # neighbours of each pixel
# here only for high gain pixels
rl = np.roll(asic_h, 1, axis=0)
rr = np.roll(asic_h, -1, axis=0)
@@ -586,10 +655,13 @@ class AgipdCorrections:
rld = np.roll(asic_h, -1, axis=1)
rru = np.roll(asic_h, 1, axis=1)
rrd = np.roll(asic_h, -1, axis=1)
- bg_h = np.nanmean(np.array([rl, rr, ru, rd, rlu, rld, rru, rrd]), axis=0)
+ bg_h = np.nanmean(
+ np.array([rl, rr, ru, rd, rlu, rld, rru, rrd]), axis=0)
# construct a graph
- graph = np.array([np.abs(asic[fidx]-bg[fidx])*asic_r[fidx]**2, asic[fidx]]).T
+ graph = np.array(
+ [np.abs(asic[fidx] - bg[fidx]) * asic_r[fidx] ** 2,
+ asic[fidx]]).T
# scale it
graph = StandardScaler().fit_transform(graph)
# perform clustering
@@ -598,14 +670,15 @@ class AgipdCorrections:
# and get labels
labels = spc.predict(graph)
- asic = im[k, i*64:(i+1)*64,j*64:(j+1)*64]
- asic_msk = snow_mask[k, i*64:(i+1)*64,j*64:(j+1)*64]
+ asic = im[k, i * 64:(i + 1) * 64, j * 64:(j + 1) * 64]
+ asic_msk = snow_mask[k, i * 64:(i + 1) * 64,
+ j * 64:(j + 1) * 64]
mim = asic[fidx]
asicb = bg_h
mimb = asicb[fidx]
mimg = asic_g[fidx]
mmsk = asic_msk[fidx]
-
+
# identify which labels are which
mn1 = np.nanmean(bg[fidx][labels == 0])
mn2 = np.nanmean(bg[fidx][labels == 1])
@@ -613,13 +686,14 @@ class AgipdCorrections:
if mn1 > mn2:
implabel = 0
- # if we've misidentified, then we will have to many snowy pixels
+ # if we've misidentified, then we will have to many
+ # snowy pixels
# happens if the ASIC is generally on a high signal level
- if np.count_nonzero([labels == implabel]) > 64*64/3:
+ if np.count_nonzero([labels == implabel]) > 64 * 64 / 3:
continue
-
+
# or a large portion of misidenfied label covers the ASIC
- if np.count_nonzero([labels != implabel]) > 0.8*64*64:
+ if np.count_nonzero([labels != implabel]) > 0.8 * 64 * 64:
continue
# set to NaN if requested
@@ -630,150 +704,171 @@ class AgipdCorrections:
mim[labels == implabel] = mimb[labels == implabel]
mimg[labels == implabel] = 0
# identify these pixels in a bad pixel mask
- mmsk[labels == implabel] = BadPixels.TRANSITION_REGION.value
-
+ mmsk[
+ labels == implabel] = BadPixels.TRANSITION_REGION.value
+
# now set back to data
asic[fidx] = mim
asic_g[fidx] = mimg
asic_msk[fidx] = mmsk
- im[k, i*64:(i+1)*64,j*64:(j+1)*64] = asic
- gain[k, i*64:(i+1)*64,j*64:(j+1)*64] = asic_g
- snow_mask[k, i*64:(i+1)*64,j*64:(j+1)*64] = asic_msk
+ im[k, i * 64:(i + 1) * 64, j * 64:(j + 1) * 64] = asic
+ gain[k, i * 64:(i + 1) * 64, j * 64:(j + 1) * 64] = asic_g
+ snow_mask[k, i * 64:(i + 1) * 64,
+ j * 64:(j + 1) * 64] = asic_msk
return im, gain, snow_mask
def correct_agipd(self, irange):
""" Correct Raw AGIOPD data for offset and relative gain effects
:param irange: list of image indices to work on, should be contiguous
-
+
Will raise an RuntimeError if `initialize()` has not been called first.
"""
if not self.initialized:
raise RuntimeError("Must call initialize() first!")
-
+
agipd_base = self.agipd_base
cidx = self.cidx
- im = np.array(self.infile[agipd_base+"image/data"][irange,...])
- trainId = np.squeeze(self.infile[agipd_base+"/image/trainId"][irange, ...])
- pulseId = np.squeeze(self.infile[agipd_base+"image/pulseId"][irange, ...])
- status = np.squeeze(self.infile[agipd_base+"/image/status"][irange, ...])
- cellid = np.squeeze(np.array(self.infile[agipd_base+"/image/cellId"][irange, ...]))
- length = np.squeeze(np.array(self.infile[agipd_base+"/image/length"][irange, ...]))
+ im = np.array(self.infile[agipd_base + "image/data"][irange, ...])
+ trainId = np.squeeze(
+ self.infile[agipd_base + "/image/trainId"][irange, ...])
+ pulseId = np.squeeze(
+ self.infile[agipd_base + "image/pulseId"][irange, ...])
+ status = np.squeeze(
+ self.infile[agipd_base + "/image/status"][irange, ...])
+ cellid = np.squeeze(
+ np.array(self.infile[agipd_base + "/image/cellId"][irange, ...]))
+ length = np.squeeze(
+ np.array(self.infile[agipd_base + "/image/length"][irange, ...]))
# split off image and gain into two separate arrays
im, ga = self.split_gain(im)
# we will work on float data from now on
im = im.astype(np.float32)
-
+
# some correction require us to maintain a copy of the raw data
raw = copy.copy(im)
-
+
# this far end of the image index range we are working on
- nidx = int(cidx+irange.size)
-
+ nidx = int(cidx + irange.size)
+
# first evaluate the gain into 0, 1, 2 --> high, medium, low
gain = np.zeros(ga.shape, np.uint8)
gain[...] = 0
- t0 = self.thresholds[...,0]
- t1 = self.thresholds[...,1]
+ t0 = self.thresholds[..., 0]
+ t1 = self.thresholds[..., 1]
# exceeding first threshold means data is medium or low gain
gain[ga > t0[cellid, ...]] = 1
# exceeding also second threshold means data is low gain
gain[ga > t1[cellid, ...]] = 2
- # check if any data has zero standard deviation, and mark this in the bad pixel maks
+ # check if any data has zero standard deviation, and mark this in
+ # the bad pixel maks
# this can be done on otherwise not corrected data.
if self.sig_zero_mask is None:
- self.sig_zero_mask = np.zeros((self.max_cells, im.shape[1], im.shape[2]), np.uint32)
+ self.sig_zero_mask = np.zeros(
+ (self.max_cells, im.shape[1], im.shape[2]), np.uint32)
for c in range(self.max_cells):
std = np.nanstd(im[cellid == c, ...], axis=0)
- self.sig_zero_mask[c, std==0] = BadPixels.DATA_STD_IS_ZERO.value
-
+ self.sig_zero_mask[
+ c, std == 0] = BadPixels.DATA_STD_IS_ZERO.value
+
# for feedback we produced histograms for the first chunk
if cidx == 0:
- H, xe, ye = np.histogram2d(im.flatten(), ga.flatten(),
- bins=self.bins_gain_vs_signal,
- range=[[4000, 8192], [4000, 8192]])
+ H, xe, ye = np.histogram2d(im.flatten(), ga.flatten(),
+ bins=self.bins_gain_vs_signal,
+ range=[[4000, 8192], [4000, 8192]])
self.hists_gain_vs_signal += H
self.signal_edges = (xe, ye)
- H, xe, ye = np.histogram2d(im.flatten(), gain.flatten(),
- bins=self.bins_dig_gain_vs_signal,
- range=[[4000, 8192], [0, 4]])
+ H, xe, ye = np.histogram2d(im.flatten(), gain.flatten(),
+ bins=self.bins_dig_gain_vs_signal,
+ range=[[4000, 8192], [0, 4]])
self.hists_dig_gain_vs_signal += H
self.dig_signal_edges = (xe, ye)
# now get the correct constants depending on cell id
- rc = self.rel_gain[cellid,...]
- offsetb = self.offset[cellid,...]
- tmask = self.mask[cellid,...]
+ rc = self.rel_gain[cellid, ...]
+ offsetb = self.offset[cellid, ...]
+ tmask = self.mask[cellid, ...]
# choose constants according to gain setting
- off = np.choose(gain, (offsetb[...,0], offsetb[...,1], offsetb[...,2]))
- rel_cor = np.choose(gain, (rc[...,0], rc[...,1], rc[...,2]))
-
+ off = np.choose(gain,
+ (offsetb[..., 0], offsetb[..., 1], offsetb[..., 2]))
+ rel_cor = np.choose(gain, (rc[..., 0], rc[..., 1], rc[..., 2]))
+
# also choose the correct bad pixel mask
- msk = np.choose(gain, (tmask[...,0], tmask[...,1], tmask[...,2]))
+ msk = np.choose(gain, (tmask[..., 0], tmask[..., 1], tmask[..., 2]))
# scale raw gain for use in the identifying snowy pixels
rgain = None
- if self.melt_snow is not False:
- rgain = ga/t0[cellid, ...]
-
+ if self.melt_snow is not False:
+ rgain = ga / t0[cellid, ...]
+
# subtract offset
im -= off
-
- # before doing relative gain correction we need to evaluate any baseline shifts
+
+ # before doing relative gain correction we need to evaluate any
+ # baseline shifts
# as they are effectively and additional offset in the data
if self.baseline_corr_using_noise or self.baseline_corr_using_hmatch:
-
+
# do this image wise, as the shift is per image
for i in range(im.shape[0]):
-
- # first correction requested may be to evaluate shift via noise peak
+
+ # first correction requested may be to evaluate shift via
+ # noise peak
if self.baseline_corr_using_noise:
- dd = self.baseline_correct_via_noise(im[i,...],
- np.nanmean(self.noise[cellid[i],...,0]),
- gain[i,...])
+ mn_noise = np.nanmean(self.noise[cellid[i], ..., 0])
+ dd = self.baseline_correct_via_noise(im[i, ...],
+ mn_noise,
+ gain[i, ...])
# if not we continue with initial data
else:
- dd = im[i,...]
-
- # if we have enough pixels in medium or low gain and correction via hist
- # matching is requested to this now
- if np.count_nonzero(gain[i,...]>0) > 1000 and self.baseline_corr_using_hmatch:
- dd2 = self.correct_baseline_via_hist(im[i,...],
- rel_cor[i,...],
- gain[i,...])
- im[i,...] = np.maximum(dd, dd2)
-
+ dd = im[i, ...]
+
+ # if we have enough pixels in medium or low gain and
+ # correction via hist matching is requested to this now
+ gcrit = np.count_nonzero(gain[i, ...] > 0) > 1000
+ if gcrit and self.baseline_corr_using_hmatch:
+ dd2 = self.correct_baseline_via_hist(im[i, ...],
+ rel_cor[i, ...],
+ gain[i, ...])
+ im[i, ...] = np.maximum(dd, dd2)
+
# finally correct diagonal effects if requested
if self.correct_asic_diag:
- im[i,...] = self.correct_baseline_via_hist_asic(im[i,...], gain[i,...])
- # if there is not enough medium or low gain data to do an evaluation, do nothing
+ ii = im[i, ...]
+ gg = gain[i, ...]
+ adim = self.correct_baseline_via_hist_asic(ii, gg)
+ im[i, ...] = adim
+ # if there is not enough medium or low gain data to do an
+ # evaluation, do nothing
else:
- im[i,...] = dd
+ im[i, ...] = dd
# now we can correct for relative gain if requested
if self.do_rel_gain:
- im *= rel_cor
-
+ im *= rel_cor
+
# try to identify snowy pixels at this point
if self.melt_snow is not False:
+ ms = self.melt_snow
im, gain, snowmask = self.melt_snowy_pixels(raw,
im,
- gain,
- rgain,
- resolution=self.melt_snow)
-
- # finally, with all corrections performed we can match ASIC borders if needed
+ gain,
+ rgain,
+ resolution=ms)
+
+ # finally, with all corrections performed we can match ASIC borders
+ # if needed
if self.match_asics:
im = self.match_asic_borders(im)
# create a bad pixel mask for the data
# we add any non-finite values to the mask
- bidx = ~np.isfinite(im)
+ bidx = ~np.isfinite(im)
im[bidx] = 0
msk[bidx] |= BadPixels.VALUE_IS_NAN.value
@@ -782,7 +877,8 @@ class AgipdCorrections:
im[bidx] = 0
msk[bidx] |= BadPixels.VALUE_OUT_OF_RANGE.value
- # include pixels with zero standard deviation in the dataset into the mask
+ # include pixels with zero standard deviation in the dataset into
+ # the mask
msk |= self.sig_zero_mask[cellid, ...]
if self.melt_snow is not False:
msk |= snowmask
@@ -791,56 +887,66 @@ class AgipdCorrections:
if cidx == 0:
copim = copy.copy(im)
copim[copim < self.median_noise] = np.nan
- H, xe, ye = np.histogram2d(np.nanmean(copim, axis=(1,2)), pulseId,
- bins=(self.bins_signal_low_range, self.max_pulses),
- range=[[-50, 1000],[0, self.max_pulses+1]])
+ bins = (self.bins_signal_low_range, self.max_pulses)
+ rnge = [[-50, 1000], [0, self.max_pulses + 1]]
+ H, xe, ye = np.histogram2d(np.nanmean(copim, axis=(1, 2)),
+ pulseId,
+ bins=bins,
+ range=rnge)
self.hists_signal_low += H
self.low_edges = (xe, ye)
- H, xe, ye = np.histogram2d(np.nanmean(copim, axis=(1,2)), pulseId,
- bins=(self.bins_signal_high_range, self.max_pulses),
- range=[[0, 200000],[0, self.max_pulses+1]])
+ bins = (self.bins_signal_high_range, self.max_pulses)
+ rnge = [[0, 200000], [0, self.max_pulses + 1]]
+ H, xe, ye = np.histogram2d(np.nanmean(copim, axis=(1, 2)),
+ pulseId,
+ bins=bins,
+ range=rnge)
self.hists_signal_high += H
self.high_edges = (xe, ye)
# now write out the data
- self.ddset[cidx:nidx,...] = im
- self.gdset[cidx:nidx,...] = gain
- self.mdset[cidx:nidx,...] = msk
+ self.ddset[cidx:nidx, ...] = im
+ self.gdset[cidx:nidx, ...] = gain
+ self.mdset[cidx:nidx, ...] = msk
- self.outfile[agipd_base+"image/cellId"][cidx:nidx] = cellid
- self.outfile[agipd_base+"image/trainId"][cidx:nidx] = trainId
- self.outfile[agipd_base+"image/pulseId"][cidx:nidx] = pulseId
- self.outfile[agipd_base+"image/status"][cidx:nidx] = status
- self.outfile[agipd_base+"image/length"][cidx:nidx] = length
- self.cidx = nidx
+ self.outfile[agipd_base + "image/cellId"][cidx:nidx] = cellid
+ self.outfile[agipd_base + "image/trainId"][cidx:nidx] = trainId
+ self.outfile[agipd_base + "image/pulseId"][cidx:nidx] = pulseId
+ self.outfile[agipd_base + "image/status"][cidx:nidx] = status
+ self.outfile[agipd_base + "image/length"][cidx:nidx] = length
+ self.cidx = nidx
def get_valid_image_idx(self):
""" Return the indices of valid data
"""
agipd_base = self.idx_base
if self.index_v == 2:
- count = np.squeeze(self.infile[agipd_base+"image/count"])
- first = np.squeeze(self.infile[agipd_base+"image/first"])
+ count = np.squeeze(self.infile[agipd_base + "image/count"])
+ first = np.squeeze(self.infile[agipd_base + "image/first"])
if np.count_nonzero(count != 0) == 0:
raise IOError("File has no valid counts")
valid = count != 0
idxtrains = np.squeeze(self.infile["/INDEX/trainId"])
medianTrain = np.nanmedian(idxtrains)
- valid &= (idxtrains > medianTrain - 1e4) & (idxtrains < medianTrain + 1e4)
+ lowok = (idxtrains > medianTrain - 1e4)
+ highok = (idxtrains < medianTrain + 1e4)
+ valid &= lowok & highok
- last_index = int(first[valid][-1]+count[valid][-1])
+ last_index = int(first[valid][-1] + count[valid][-1])
first_index = int(first[valid][0])
elif self.index_v == 1:
- status = np.squeeze(self.infile[agipd_base+"image/status"])
+ status = np.squeeze(self.infile[agipd_base + "image/status"])
if np.count_nonzero(status != 0) == 0:
- raise IOError("File {} has no valid counts".format(infile))
- last = np.squeeze(self.infile[agipd_base+"image/last"])
+ raise IOError("File {} has no valid counts".format(
+ self.infile))
+ last = np.squeeze(self.infile[agipd_base + "image/last"])
valid = status != 0
last_index = int(last[status != 0][-1])
first_index = int(last[status != 0][0])
else:
- raise AttributeError("Not a known raw format version: {}".format(self.index_v))
+ raise AttributeError(
+ "Not a known raw format version: {}".format(self.index_v))
self.valid = valid
self.first_index = first_index
@@ -854,20 +960,26 @@ class AgipdCorrections:
last_index = self.last_index
max_cells = self.max_cells
agipd_base = self.agipd_base
- allcells = np.squeeze(np.array(self.infile[agipd_base+"image/cellId"][first_index:last_index, ...]))
- single_image = np.array(np.array(self.infile[agipd_base+"image/data"][first_index, ...]))
+ allcells = self.infile[agipd_base + "image/cellId"]
+ allcells = np.squeeze(allcells[first_index:last_index, ...])
+ single_image = self.infile[agipd_base + "image/data"][first_index, ...]
+ single_image = np.array(single_image)
+
can_calibrate = allcells < max_cells
if np.count_nonzero(can_calibrate) == 0:
return
allcells = allcells[can_calibrate]
firange = np.arange(first_index, last_index)
firange = firange[can_calibrate]
- self.oshape = (firange.size, single_image.shape[1], single_image.shape[2])
+ self.oshape = (firange.size,
+ single_image.shape[1],
+ single_image.shape[2])
self.firange = firange
self.single_image = single_image
def copy_and_sanitize_non_cal_data(self):
- """ Copy and sanitize data in `infile` that is not touched by `correctAGIPD`
+ """ Copy and sanitize data in `infile` that is not touched by
+ `correctAGIPD`
"""
agipd_base = self.agipd_base
idx_base = self.idx_base
@@ -876,18 +988,20 @@ class AgipdCorrections:
valid = self.valid
idxtrains = self.idxtrains
firange = self.firange
- alltrains = np.squeeze(self.infile[agipd_base+"image/trainId"][first_index:last_index,...])
-
+ alltrains = self.infile[agipd_base + "image/trainId"]
+ alltrains = np.squeeze(alltrains[first_index:last_index, ...])
+
# these are touched in the correct function, do not copy them here
- dont_copy = ["data", "cellId", "trainId", "pulseId", "status", "length"]
- dont_copy = [agipd_base+"image/{}".format(do)
- for do in dont_copy]
-
+ dont_copy = ["data", "cellId", "trainId", "pulseId", "status",
+ "length"]
+ dont_copy = [agipd_base + "image/{}".format(do)
+ for do in dont_copy]
+
# don't copy these as we may need to adjust if we filter trains
- dont_copy += [idx_base+"{}/first".format(do)
- for do in ["image",]]
- dont_copy += [idx_base+"{}/count".format(do)
- for do in ["image",]]
+ dont_copy += [idx_base + "{}/first".format(do)
+ for do in ["image", ]]
+ dont_copy += [idx_base + "{}/count".format(do)
+ for do in ["image", ]]
# a visitor to copy everything else
def visitor(k, item):
@@ -901,19 +1015,23 @@ class AgipdCorrections:
self.infile.copy(k, self.outfile[group])
self.infile.visititems(visitor)
-
+
# sanitize indices
for do in ["image", ]:
- existing = np.squeeze(self.infile[idx_base+"{}/first".format(do)])
- updated = existing-np.min(existing[valid])
- self.outfile[idx_base+"{}/first".format(do)] = updated
-
- existing = np.squeeze(self.infile[idx_base+"{}/count".format(do)])
- new_counts, _ = np.histogram(alltrains[firange], bins=np.count_nonzero(valid)+1,
- range=(np.min(idxtrains[valid]), np.max(idxtrains[valid]+1)))
+ existing = np.squeeze(
+ self.infile[idx_base + "{}/first".format(do)])
+ updated = existing - np.min(existing[valid])
+ self.outfile[idx_base + "{}/first".format(do)] = updated
+
+ existing = np.squeeze(
+ self.infile[idx_base + "{}/count".format(do)])
+ new_counts, _ = np.histogram(alltrains[firange],
+ bins=np.count_nonzero(valid) + 1,
+ range=(np.min(idxtrains[valid]),
+ np.max(idxtrains[valid] + 1)))
updated = existing
updated[valid] = new_counts[:-1]
- self.outfile[idx_base+"{}/count".format(do)] = updated
+ self.outfile[idx_base + "{}/count".format(do)] = updated
def create_output_datasets(self):
""" Initialize output data sets
@@ -922,279 +1040,335 @@ class AgipdCorrections:
chunksize = self.chunk_size_idim
firange = self.firange
oshape = self.oshape
- self.ddset = self.outfile.create_dataset(agipd_base+"image/data",
- oshape, chunks=(chunksize, oshape[1], oshape[2]), dtype=np.float32)
- self.gdset = self.outfile.create_dataset(agipd_base+"image/gain",
- oshape, chunks=(chunksize, oshape[1], oshape[2]), dtype=np.uint8,
- compression="gzip", compression_opts=1, shuffle=True)
- self.mdset = self.outfile.create_dataset(agipd_base+"image/mask",
- oshape, chunks=(chunksize, oshape[1], oshape[2]), dtype=np.uint32,
- compression="gzip", compression_opts=1, shuffle=True)
-
- self.outfile[agipd_base+"image/cellId"] = np.zeros(firange.size, np.uint16)
- self.outfile[agipd_base+"image/trainId"] = np.zeros(firange.size, np.uint64)
- self.outfile[agipd_base+"image/pulseId"] = np.zeros(firange.size, np.uint64)
- self.outfile[agipd_base+"image/status"] = np.zeros(firange.size, np.uint16)
- self.outfile[agipd_base+"image/length"] = np.zeros(firange.size, np.uint32)
+ chunks = (chunksize, oshape[1], oshape[2])
+ self.ddset = self.outfile.create_dataset(agipd_base + "image/data",
+ oshape, chunks=chunks,
+ dtype=np.float32)
+ self.gdset = self.outfile.create_dataset(agipd_base + "image/gain",
+ oshape, chunks=chunks,
+ dtype=np.uint8,
+ compression="gzip",
+ compression_opts=1,
+ shuffle=True)
+ self.mdset = self.outfile.create_dataset(agipd_base + "image/mask",
+ oshape, chunks=chunks,
+ dtype=np.uint32,
+ compression="gzip",
+ compression_opts=1,
+ shuffle=True)
+
+ fsz = firange.size
+ self.outfile[agipd_base + "image/cellId"] = np.zeros(fsz, np.uint16)
+ self.outfile[agipd_base + "image/trainId"] = np.zeros(fsz, np.uint64)
+ self.outfile[agipd_base + "image/pulseId"] = np.zeros(fsz, np.uint64)
+ self.outfile[agipd_base + "image/status"] = np.zeros(fsz, np.uint16)
+ self.outfile[agipd_base + "image/length"] = np.zeros(fsz, np.uint32)
self.outfile.flush()
-
+
def get_histograms(self):
""" Return preview histograms computed from the first chunk
"""
- return ((self.hists_signal_low, self.hists_signal_high, self.hists_gain_vs_signal, self.hists_dig_gain_vs_signal),
- (self.low_edges, self.high_edges, self.signal_edges, self.dig_signal_edges))
-
+ return ((self.hists_signal_low, self.hists_signal_high,
+ self.hists_gain_vs_signal, self.hists_dig_gain_vs_signal),
+ (self.low_edges, self.high_edges, self.signal_edges,
+ self.dig_signal_edges))
+
def initialize_from_db(self, dbparms, qm, only_dark=False):
""" Initialize calibration constants from the calibration database
-
- :param dbparms: a tuple containing relevant database parameters, can be either:
- * cal_db_interface, creation_time, max_cells_db, bias_voltage, photon_energy
- in which case the db timeout is set to 300 seconds, the the cells to query
- dark image derived constants from the database is set to the global value
-
- * cal_db_interface, creation_time, max_cells_db, bias_voltage, photon_energy, max_cells_db_dark
- here the number of memory cells to query dark derived data differs from the global value, the
- the db timeout is set to 300 seconds
-
- * cal_db_interface, creation_time, max_cells_db, bias_voltage, photon_energy,
- max_cells_db_dark, timeout_cal_db
- addiitonally a timeout is given
-
- :param qm: quadrant and module of the constants to load in Q1M1 notation
- :param only_dark: load only dark image derived constants. This implies that a `calfile` is
- used to load the remaining constants. Useful to reduce DB traffic and interactions
- for non-frequently changing constants, i.e. such which are not usually updated during
- a beamtime.
-
- The `cal_db_interface` parameter in the `dbparms` tuple may be in one of the following notations:
- * tcp://host:port to directly identify the host and port to connect to
- * tcp://host:port_low#port_high to specify a port range from which a random port will be picked.
-
+
+ :param dbparms: a tuple containing relevant database parameters,
+ can be either:
+ * cal_db_interface, creation_time, max_cells_db, bias_voltage,
+ photon_energy in which case the db timeout is set to 300 seconds,
+ the cells to query dark image derived constants from the
+ database is set to the global value
+
+ * cal_db_interface, creation_time, max_cells_db, bias_voltage,
+ photon_energy, max_cells_db_dark here the number of memory
+ cells to query dark derived data differs from the global
+ value, the the db timeout is set to 300 seconds
+
+ * cal_db_interface, creation_time, max_cells_db, bias_voltage,
+ photon_energy, max_cells_db_dark, timeout_cal_db
+ additionally a timeout is given
+
+ :param qm: quadrant and module of the constants to load in Q1M1
+ notation
+ :param only_dark: load only dark image derived constants. This
+ implies that a `calfile` is used to load the remaining
+ constants. Useful to reduce DB traffic and interactions
+ for non-frequently changing constants, i.e. such which are
+ not usually updated during a beamtime.
+
+ The `cal_db_interface` parameter in the `dbparms` tuple may be in
+ one of the following notations:
+ * tcp://host:port to directly identify the host and port to
+ connect to
+ * tcp://host:port_low#port_high to specify a port range from
+ which a random port will be picked.
+
The latter notation allows for load-balancing.
-
- This routine loads the following constants as given in `iCalibrationDB`:
-
+
+ This routine loads the following constants as given in
+ `iCalibrationDB`:
+
Dark Image Derived
------------------
-
+
* Constants.AGIPD.Offset
* Constants.AGIPD.Noise
* Constants.AGIPD.BadPixelsDark
* Constants.AGIPD.ThresholdsDark
-
+
Pulse Capacitor Derived
-----------------------
-
- * Constants.AGIPD.SlopesPC - older data, whith the parameter dimension first is correctly handled
-
+
+ * Constants.AGIPD.SlopesPC - older data, whith the parameter
+ dimension first is correctly handled
+
Flat-Field Derived
-
+
* Constants.AGIPD.SlopesFF
-
- For CI derived gain, a mean multiplication factor of 4.48 compared to medium gain is used, as no
- reliable CI data for all memory cells exists for all memory cells of the current AGIPD instances.
-
- Relative gain is derived both from pulse capacitor as well as flat field data:
-
- * from the pulse capacitor data we get the relative slopes of a given pixels memory cells
- with respect to all memory cells of that pixel:
-
+
+ For CI derived gain, a mean multiplication factor of 4.48 compared
+ to medium gain is used, as no reliable CI data for all memory cells
+ exists for all memory cells of the current AGIPD instances.
+
+ Relative gain is derived both from pulse capacitor as well as flat
+ field data:
+
+ * from the pulse capacitor data we get the relative slopes of a
+ given pixels memory cells with respect to all memory cells of that
+ pixel:
+
rpch = m_h / median(m_h)
-
+
where m_h is the high gain slope m of each memory cell of the pixel.
-
- * we also derive the factor between high and medium gain in a similar way and scale it to be relative
- to the pixels memory cells:
-
+
+ * we also derive the factor between high and medium gain in a
+ similar way and scale it to be relative to the pixels memory cells:
+
fpc = m_m/m_h
rfpc = fpc/ median(fpc)
-
- where m_m is the medium gain slope of all memory cells of a given pixel and m_h is the high gain
- slope as above
-
- * finally, we get the relative X-ray gain of all memory cells for a given pixel from flat field
- data:
-
+
+ where m_m is the medium gain slope of all memory cells of a given
+ pixel and m_h is the high gain slope as above
+
+ * finally, we get the relative X-ray gain of all memory cells for a
+ given pixel from flat field data:
+
ff = median(m_ff)
ff /= median(ff)
-
- where m_ff is the flat field derived (high gain) slope of all memory cells of a given pixel.
- The pixel values are then scaled to the complete module. Note that the first 32 memory cells
- are known to exhibit differing behaviour and are skipped if possible.
-
- With this data the relative gain for the three gain stages evaluates to:
-
+
+ where m_ff is the flat field derived (high gain) slope of all
+ memory cells of a given pixel. The pixel values are then scaled to
+ the complete module. Note that the first 32 memory cells are known
+ to exhibit differing behaviour and are skipped if possible.
+
+ With this data the relative gain for the three gain stages evaluates
+ to:
+
high gain = ff * rpch
medium gain = ff * rfpc
low gain = medium gain / 4.48
-
+
"""
if len(dbparms) == 5:
- cal_db_interface, creation_time, max_cells_db, bias_voltage, photon_energy = dbparms
+ (cal_db_interface, creation_time, max_cells_db,
+ bias_voltage, photon_energy) = dbparms
max_cells_db_dark = max_cells_db
timeout_cal_db = 300000
-
+
elif len(dbparms) == 6:
(cal_db_interface, creation_time, max_cells_db,
bias_voltage, photon_energy, max_cells_db_dark) = dbparms
-
+
if max_cells_db_dark is None:
max_cells_db_dark = max_cells_db
timeout_cal_db = 300000
else:
(cal_db_interface, creation_time, max_cells_db,
- bias_voltage, photon_energy, max_cells_db_dark, timeout_cal_db) = dbparms
-
+ bias_voltage, photon_energy, max_cells_db_dark,
+ timeout_cal_db) = dbparms
+
if max_cells_db_dark is None:
max_cells_db_dark = max_cells_db
-
+
if "#" in cal_db_interface:
prot, serv, ran = cal_db_interface.split(":")
r1, r2 = ran.split("#")
- cal_db_interface = ":".join([prot, serv, str(np.random.randint(int(r1), int(r2)))])
-
+ cal_db_interface = ":".join(
+ [prot, serv, str(np.random.randint(int(r1), int(r2)))])
+
offset = get_constant_from_db(getattr(Detectors.AGIPD1M1, qm),
Constants.AGIPD.Offset(),
- Conditions.Dark.AGIPD(memory_cells=max_cells_db_dark,
- bias_voltage=bias_voltage),
- np.zeros((128,512,max_cells_db,3)),
+ Conditions.Dark.AGIPD(
+ memory_cells=max_cells_db_dark,
+ bias_voltage=bias_voltage),
+ np.zeros((128, 512, max_cells_db, 3)),
cal_db_interface,
creation_time=creation_time,
timeout=timeout_cal_db)
-
+
noise = get_constant_from_db(getattr(Detectors.AGIPD1M1, qm),
Constants.AGIPD.Noise(),
- Conditions.Dark.AGIPD(memory_cells=max_cells_db_dark,
- bias_voltage=bias_voltage),
- np.zeros((128,512,max_cells_db,3)),
+ Conditions.Dark.AGIPD(
+ memory_cells=max_cells_db_dark,
+ bias_voltage=bias_voltage),
+ np.zeros((128, 512, max_cells_db, 3)),
cal_db_interface,
creation_time=creation_time,
timeout=timeout_cal_db)
-
+
bpixels = get_constant_from_db(getattr(Detectors.AGIPD1M1, qm),
Constants.AGIPD.BadPixelsDark(),
- Conditions.Dark.AGIPD(memory_cells=max_cells_db_dark,
- bias_voltage=bias_voltage),
- np.zeros((128,512,max_cells_db,3)),
+ Conditions.Dark.AGIPD(
+ memory_cells=max_cells_db_dark,
+ bias_voltage=bias_voltage),
+ np.zeros((128, 512, max_cells_db, 3)),
cal_db_interface,
creation_time=creation_time,
- timeout=timeout_cal_db).astype(np.uint32)
-
+ timeout=timeout_cal_db).astype(
+ np.uint32)
+
thresholds = get_constant_from_db(getattr(Detectors.AGIPD1M1, qm),
Constants.AGIPD.ThresholdsDark(),
- Conditions.Dark.AGIPD(memory_cells=max_cells_db_dark,
- bias_voltage=bias_voltage),
- np.ones((128,512,max_cells_db,2)),
+ Conditions.Dark.AGIPD(
+ memory_cells=max_cells_db_dark,
+ bias_voltage=bias_voltage),
+ np.ones((128, 512, max_cells_db, 2)),
cal_db_interface,
creation_time=creation_time,
timeout=timeout_cal_db)
-
+
# we have all we need for dark image derived constants
# initialize and return
if only_dark:
self.initialize(offset, None, bpixels, noise, thresholds, None)
- return
+ return
- # load also non-dark constants from the database
+ # load also non-dark constants from the database
slopesPC = get_constant_from_db(getattr(Detectors.AGIPD1M1, qm),
Constants.AGIPD.SlopesPC(),
- Conditions.Dark.AGIPD(memory_cells=max_cells_db,
- bias_voltage=bias_voltage),
- np.ones((128,512,max_cells_db,10)),
+ Conditions.Dark.AGIPD(
+ memory_cells=max_cells_db,
+ bias_voltage=bias_voltage),
+ np.ones((128, 512, max_cells_db, 10)),
cal_db_interface,
creation_time=creation_time,
- timeout=timeout_cal_db).astype(np.float32)
-
+ timeout=timeout_cal_db)
+ slopesPC = slopesPC.astype(np.float32)
+
# this will handle some historical data in a slighly different format
- if slopesPC.shape[0] == 10: # constant dimension injected first
+ if slopesPC.shape[0] == 10: # constant dimension injected first
slopesPC = np.moveaxis(slopesPC, 0, 3)
slopesPC = np.moveaxis(slopesPC, 0, 2)
-
- slopesFF = np.squeeze(get_constant_from_db(getattr(Detectors.AGIPD1M1, qm),
- Constants.AGIPD.SlopesFF(),
- Conditions.Illuminated.AGIPD(max_cells_db, bias_voltage, photon_energy,
- pixels_x=512, pixels_y=128,
- beam_energy=None),
- np.ones((128,512,max_cells_db,2)),
- cal_db_interface,
- creation_time=creation_time,
- timeout=timeout_cal_db))
-
+
+ slopesFF = np.squeeze(
+ get_constant_from_db(getattr(Detectors.AGIPD1M1, qm),
+ Constants.AGIPD.SlopesFF(),
+ Conditions.Illuminated.AGIPD(max_cells_db,
+ bias_voltage,
+ photon_energy,
+ pixels_x=512,
+ pixels_y=128,
+ beam_energy=None), # noqa
+ np.ones((128, 512, max_cells_db, 2)),
+ cal_db_interface,
+ creation_time=creation_time,
+ timeout=timeout_cal_db))
+
# calculate relative gain from the constants
# we default to a value of one, so basically a unity transform
rel_gain = np.ones((128, 512, self.max_cells, 3), np.float32)
-
+
# high gain slope from pulse capacitor data
- pc_high_m = slopesPC[...,:self.max_cells,0]
-
+ pc_high_m = slopesPC[..., :self.max_cells, 0]
+
# factor of high vs. medium gain for each pixel as derived from PC data
- fac_high_med = slopesPC[...,:self.max_cells,3] / slopesPC[...,:self.max_cells,0]
-
- # if the flat field data contains a field for X-ray offset we can make used of it
+ fac_high_med = slopesPC[..., :self.max_cells, 3] / pc_high_m
+
+ # if the flat field data contains a field for X-ray offset we can
+ # make used of it
xrayOffset = 0
if len(slopesFF.shape) == 4:
- xrayOffset = slopesFF[...,:self.max_cells,1]
- slopesFF = slopesFF[...,0]
-
- # first 32 cells are known to behave differently so if we can avoid them
+ xrayOffset = slopesFF[..., :self.max_cells, 1]
+ slopesFF = slopesFF[..., 0]
+
+ # first 32 cells are known to behave differently so if we can avoid
+ # them
# when calculating the mean X-ray derived gain slope for each pixel
if slopesFF.shape[2] > 32:
- xraygain = np.nanmedian(slopesFF[...,32:min(96,self.max_cells)], axis=2)
+ xraygain = np.nanmedian(slopesFF[..., 32:min(96, self.max_cells)],
+ axis=2)
else:
- xraygain = np.nanmedian(slopesFF[...,:min(96,self.max_cells)], axis=2)
-
- # relative X-ray gain is then each pixels gain, by the median gain of all pixels
+ xraygain = np.nanmedian(slopesFF[..., :min(96, self.max_cells)],
+ axis=2)
+
+ # relative X-ray gain is then each pixels gain, by the median gain
+ # of all pixels
xraygain /= np.nanmedian(xraygain)
-
- # scale the relative PC gain of each memory cell in a pixel to the median value of all cells
- pcrel = pc_high_m/np.nanmedian(pc_high_m, axis=2)[...,None]
-
- # and get high gain relative gain by multiplying with x-ray derived gain
- rel_gain[..., 0] = pcrel*xraygain[...,None]
-
+
+ # scale the relative PC gain of each memory cell in a pixel to the
+ # median value of all cells
+ pcrel = pc_high_m / np.nanmedian(pc_high_m, axis=2)[..., None]
+
+ # and get high gain relative gain by multiplying with x-ray derived
+ # gain
+ rel_gain[..., 0] = pcrel * xraygain[..., None]
+
# same scaling for high-to medium gain factor
- mfac = fac_high_med/np.nanmedian(fac_high_med, axis=2)[...,None]*np.nanmedian(fac_high_med)
-
+ medfhm = np.nanmedian(fac_high_med, axis=2)[..., None]
+ mfac = fac_high_med / medfhm * np.nanmedian(fac_high_med)
+
# and scale with X-ray derived gain
- rel_gain[..., 1] = xraygain[...,None] * mfac
-
+ rel_gain[..., 1] = xraygain[..., None] * mfac
+
# finally low gain with a constant factor
rel_gain[..., 2] = rel_gain[..., 1] * 0.223
-
+
# for base offsets we pass a tuple containing the relevant information
- base_offset = [slopesPC[...,:self.max_cells,1], slopesPC[...,:self.max_cells,4], xrayOffset]
-
+ base_offset = [slopesPC[..., :self.max_cells, 1],
+ slopesPC[..., :self.max_cells, 4], xrayOffset]
+
# add additonal bad pixel information
bppc = get_constant_from_db(getattr(Detectors.AGIPD1M1, qm),
Constants.AGIPD.BadPixelsPC(),
- Conditions.Dark.AGIPD(memory_cells=max_cells_db,
- bias_voltage=bias_voltage),
- np.zeros((max_cells_db, 128,512)),
+ Conditions.Dark.AGIPD(
+ memory_cells=max_cells_db,
+ bias_voltage=bias_voltage),
+ np.zeros((max_cells_db, 128, 512)),
cal_db_interface,
creation_time=creation_time,
timeout=timeout_cal_db)
- bpixels |= np.moveaxis(bppc.astype(np.uint32), 0, 2)[...,:bpixels.shape[2], None]
-
+ bppc = np.moveaxis(bppc.astype(np.uint32), 0, 2)
+ bpixels |= bppc[..., :bpixels.shape[2], None]
+
bpff = get_constant_from_db(getattr(Detectors.AGIPD1M1, qm),
Constants.AGIPD.BadPixelsFF(),
- Conditions.Illuminated.AGIPD(max_cells_db, bias_voltage, photon_energy,
- pixels_x=512, pixels_y=128,
- beam_energy=None),
- np.zeros((128,512,max_cells_db)),
+ Conditions.Illuminated.AGIPD(max_cells_db,
+ bias_voltage,
+ photon_energy,
+ pixels_x=512,
+ pixels_y=128,
+ beam_energy=None), # noqa
+ np.zeros((128, 512, max_cells_db)),
cal_db_interface,
creation_time=creation_time,
timeout=timeout_cal_db)
-
- bpixels |= bpff.astype(np.uint32)[...,:bpixels.shape[2], None]
- self.initialize(offset, rel_gain, bpixels, noise, thresholds, base_offset)
-
+ bpixels |= bpff.astype(np.uint32)[..., :bpixels.shape[2], None]
+
+ self.initialize(offset, rel_gain, bpixels, noise, thresholds,
+ base_offset)
+
def initialize_from_file(self, filename, qm, with_dark=True):
""" Initialize calibration constants from a calibration file
-
- :param filename: path to a file containing the calibration constants. It is
+
+ :param filename: path to a file containing the calibration
+ constants. It is
expected to have the following structure:
-
+
/{qm}/BadPixelsFF/data
/{qm}/BadPixelsPC/data
/{qm}/Offset/data
@@ -1203,97 +1377,106 @@ class AgipdCorrections:
/{qm}/BadPixelsDark/data
/{qm}/SlopesFF/data
/{qm}/SlopesPC/data
-
+
where qm is the `qm` parameter.
-
- :param qm: quadrant and module of the constants to load in Q1M1 notation
- :param with_dark: also load dark image derived constants from the file. This will overwrite
- any constants previously loaded from the calibration database.
-
- For CI derived gain, a mean multiplication factor of 4.48 compared to medium gain is used, as no
- reliable CI data for all memory cells exists for all memory cells of the current AGIPD instances.
-
- Relative gain is derived both from pulse capacitor as well as flat field data:
-
- * from the pulse capacitor data we get the relative slopes of a given pixels memory cells
- with respect to all memory cells of that pixel:
-
+
+ :param qm: quadrant and module of the constants to load in Q1M1
+ notation
+ :param with_dark: also load dark image derived constants from the
+ file. This will overwrite any constants previously loaded from
+ the calibration database.
+
+ For CI derived gain, a mean multiplication factor of 4.48 compared
+ to medium gain is used, as no reliable CI data for all memory cells
+ exists for all memory cells of the current AGIPD instances.
+
+ Relative gain is derived both from pulse capacitor as well as flat
+ field data:
+
+ * from the pulse capacitor data we get the relative slopes of a given
+ pixels memory cells with respect to all memory cells of that pixel:
+
rpch = m_h / median(m_h)
-
+
where m_h is the high gain slope m of each memory cell of the pixel.
-
- * we also derive the factor between high and medium gain in a similar way and scale it to be relative
- to the pixels memory cells:
-
+
+ * we also derive the factor between high and medium gain in a similar
+ way and scale it to be relative to the pixels memory cells:
+
fpc = m_m/m_h
rfpc = fpc/ median(fpc)
-
- where m_m is the medium gain slope of all memory cells of a given pixel and m_h is the high gain
- slope as above
-
- * finally, we get the relative X-ray gain of all memory cells for a given pixel from flat field
- data:
-
+
+ where m_m is the medium gain slope of all memory cells of a given
+ pixel and m_h is the high gain slope as above
+
+ * finally, we get the relative X-ray gain of all memory cells for a
+ given pixel from flat field data:
+
ff = median(m_ff)
ff /= median(ff)
-
- where m_ff is the flat field derived (high gain) slope of all memory cells of a given pixel.
- The pixel values are then scaled to the complete module. Note that the first 32 memory cells
- are known to exhibit differing behaviour and are skipped if possible.
-
- With this data the relative gain for the three gain stages evaluates to:
-
+
+ where m_ff is the flat field derived (high gain) slope of all memory
+ cells of a given pixel. The pixel values are then scaled to the
+ complete module. Note that the first 32 memory cells are known to
+ exhibit differing behaviour and are skipped if possible.
+
+ With this data the relative gain for the three gain stages evaluates
+ to:
+
high gain = ff * rpch
medium gain = ff * rfpc
low gain = medium gain / 4.48
-
+
"""
offsets = None
- rel_gains = None
bpixels = None
noises = None
thresholds = None
with h5py.File(filename, "r") as calfile:
-
+
bpixels = calfile["{}/{}/data".format(qm, "BadPixelsFF")][()]
- bpixels |= np.moveaxis(calfile["{}/{}/data".format(qm, "BadPixelsPC")][()], 0, 2)
- bpixels = bpixels[...,None] # without gain dimension up to here
+ bpixels |= np.moveaxis(
+ calfile["{}/{}/data".format(qm, "BadPixelsPC")][()], 0, 2)
+ bpixels = bpixels[..., None] # without gain dimension up to here
if with_dark:
offsets = calfile["{}/{}/data".format(qm, "Offset")][()]
noises = calfile["{}/{}/data".format(qm, "Noise")][()]
- thresholds = calfile["{}/{}/data".format(qm, "ThresholdsDark")][()]
- bpixels |= calfile["{}/{}/data".format(qm, "BadPixelsDark")][()]
-
+ thresholds = calfile["{}/{}/data".format(qm, "ThresholdsDark")][()] # noqa
+ bpixels |= calfile["{}/{}/data".format(qm, "BadPixelsDark")][
+ ()]
+
slopesFF = calfile["{}/{}/data".format(qm, "SlopesFF")][()]
slopesPC = calfile["{}/{}/data".format(qm, "SlopesPC")][()]
- if slopesPC.shape[0] == 10: # constant dimension injected first
+ if slopesPC.shape[0] == 10: # constant dimension injected first
slopesPC = np.moveaxis(slopesPC, 0, 3)
slopesPC = np.moveaxis(slopesPC, 0, 2)
-
+
rel_gain = np.ones((128, 512, self.max_cells, 3), np.float32)
- pc_high_m = slopesPC[...,:self.max_cells,0]
- pc_med_m = slopesPC[...,:self.max_cells,3]
- fac_high_med = slopesPC[...,:self.max_cells,3] / slopesPC[...,:self.max_cells,0]
-
+ pc_high_m = slopesPC[..., :self.max_cells, 0]
+ pc_med_m = slopesPC[..., :self.max_cells, 3]
+ fac_high_med = pc_med_m / pc_high_m
+
xrayOffset = 0
if len(slopesFF.shape) == 4:
- xrayOffset = slopesFF[...,:self.max_cells,1]
- slopesFF = slopesFF[...,0]
+ xrayOffset = slopesFF[..., :self.max_cells, 1]
+ slopesFF = slopesFF[..., 0]
if slopesFF.shape[2] > 32:
- xraygain = np.nanmedian(slopesFF[...,32:], axis=2)
+ xraygain = np.nanmedian(slopesFF[..., 32:], axis=2)
else:
- xraygain = np.nanmedian(slopesFF[...,:min(96,self.max_cells)], axis=2)
-
+ xraygain = np.nanmedian(
+ slopesFF[..., :min(96, self.max_cells)], axis=2)
+
xraygain /= np.nanmedian(xraygain)
- pcrel = pc_high_m/np.nanmedian(pc_high_m, axis=2)[...,None]
-
- rel_gain[..., 0] = pcrel*xraygain[...,None]
- mfac = fac_high_med/np.nanmedian(fac_high_med, axis=2)[...,None]*np.nanmedian(fac_high_med)
- rel_gain[..., 1] = xraygain[...,None] * mfac
+ pcrel = pc_high_m / np.nanmedian(pc_high_m, axis=2)[..., None]
+
+ rel_gain[..., 0] = pcrel * xraygain[..., None]
+ mfac = fac_high_med / np.nanmedian(fac_high_med, axis=2)[
+ ..., None] * np.nanmedian(fac_high_med)
+ rel_gain[..., 1] = xraygain[..., None] * mfac
rel_gain[..., 2] = rel_gain[..., 1] * 0.223
-
- base_offset = [slopesPC[...,:self.max_cells,1], slopesPC[...,:self.max_cells,4], xrayOffset]
-
- self.initialize(offsets, rel_gain, bpixels, noises, thresholds, base_offset)
-
\ No newline at end of file
+ base_offset = [slopesPC[..., :self.max_cells, 1],
+ slopesPC[..., :self.max_cells, 4], xrayOffset]
+
+ self.initialize(offsets, rel_gain, bpixels, noises, thresholds,
+ base_offset)
--
GitLab