diff --git a/cal_tools/cal_tools/agipdlib.py b/cal_tools/cal_tools/agipdlib.py
index 5d0263b11bc5fea27327891d607b9cdf600fc3de..919025af8dfee77260fc494f6434869b21512d1e 100644
--- a/cal_tools/cal_tools/agipdlib.py
+++ b/cal_tools/cal_tools/agipdlib.py
@@ -2,23 +2,23 @@ import copy
from enum import Enum
import h5py
+from iCalibrationDB import Constants, Conditions, Detectors
import numpy as np
from scipy.signal import cwt, ricker
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler
from cal_tools.enums import BadPixels
-from cal_tools.tools import get_constant_from_db
-from iCalibrationDB import Constants, Conditions, Detectors
+from cal_tools.tools import get_constant_from_db, get_constant_from_db_and_time
def get_num_cells(fname, loc, module):
with h5py.File(fname, "r") as f:
-
- cells = f["INSTRUMENT/{}/DET/{}CH0:xtdf/image/cellId".format(loc, module)][:352]
+ cells = \
+ f["INSTRUMENT/{}/DET/{}CH0:xtdf/image/cellId".format(loc, module)][()]
maxcell = np.max(cells)
options = [4, 32, 64, 76, 128, 176]
- dists = [abs(o-maxcell) for o in options]
+ dists = [abs(o - maxcell) for o in options]
return options[np.argmin(dists)]
@@ -65,7 +65,8 @@ class AgipdCorrections:
h5_index_path="INDEX/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/",
do_rel_gain=True, chunk_size_idim=512, il_mode=False,
cal_det_instance="AGIPD1M1", karabo_data_mode=False,
- force_hg_if_below=None, force_mg_if_below=None, mask_noisy_adc=False):
+ force_hg_if_below=None, force_mg_if_below=None,
+ mask_noisy_adc=False, adjust_mg_baseline=False):
"""
Initialize an AgipdCorrections Class
@@ -134,16 +135,18 @@ class AgipdCorrections:
self.mask_noisy_adc = mask_noisy_adc
self.adc_mask = None
self.gain_stats = [0, 0, 0]
+ self.adjust_mg_baseline = adjust_mg_baseline
+ self.mg_bl_adjust = 0
def get_iteration_range(self):
"""Returns a range expression over which to iterate in chunks
"""
splitlength = self.firange.size // self.chunksize
splits = np.arange(0, self.firange.size, self.chunksize)
- return [ a for a in np.array_split(self.firange, splits) if a.size > 0]
+ return [a for a in np.array_split(self.firange, splits) if a.size > 0]
def initialize(self, offset, rel_gain, mask, noise, thresholds,
- base_offset, swap_axis=True):
+ base_offset, slopesPC=None, 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
@@ -183,8 +186,13 @@ class AgipdCorrections:
mvd = np.moveaxis(np.moveaxis(noise, 2, 0), 2, 1)
noise = np.ascontiguousarray(mvd)
if base_offset is not None:
+ fallback_bo = np.nanmedian(np.array(
+ [b for b in base_offset if isinstance(b, np.ndarray)]),
+ axis=0)
for i in range(len(base_offset)):
if base_offset[i] is not None:
+ if not isinstance(base_offset[i], np.ndarray):
+ base_offset[i] = fallback_bo
mvd = np.moveaxis(base_offset[i], 2, 0)
mvd = np.moveaxis(mvd, 2, 1)
base_offset[i] = np.ascontiguousarray(mvd)
@@ -221,7 +229,7 @@ class AgipdCorrections:
for i in range(8):
for j in range(2):
asico = self.offset[:, i * 64:(i + 1) * 64,
- j * 64:(j + 1) * 64, 0]
+ j * 64:(j + 1) * 64, 0]
asic_m = np.nanmedian(asico, axis=(1, 2))
global_m = np.nanmedian(self.offset[:, ..., 0],
axis=(1, 2))
@@ -241,6 +249,27 @@ class AgipdCorrections:
self.offset[..., 1] += doff
+ if self.adjust_mg_baseline and slopesPC is not None:
+ x = np.linspace(0, 1000, 1000)
+ m_h = np.moveaxis(
+ np.moveaxis(slopesPC[..., :self.max_cells, 0], 0, 2), 0, 1)
+ b_h = np.moveaxis(
+ np.moveaxis(slopesPC[..., :self.max_cells, 1], 0, 2), 0, 1)
+ m_m = np.moveaxis(
+ np.moveaxis(slopesPC[..., :self.max_cells, 3], 0, 2), 0, 1)
+ b_m = np.moveaxis(
+ np.moveaxis(slopesPC[..., :self.max_cells, 4], 0, 2), 0, 1)
+ b_ho = self.offset[..., 0]
+ b_mo = self.offset[..., 1]
+ mz = -(b_m - b_mo) / m_m
+ cx = np.zeros(mz.shape)
+ for i in range(cx.shape[0]):
+ for j in range(cx.shape[1]):
+ for k in range(cx.shape[2]):
+ cx[i, j, k] = x[np.argmin(np.abs(mz[i, j, k] - x))]
+ doz = (m_h * cx + b_h - m_m * cx - b_m)
+ self.mg_bl_adjust = doz
+
# 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:
@@ -272,7 +301,7 @@ class AgipdCorrections:
""" Split gain information off AGIPD Data
"""
return d[:, 0, ...], d[:, 1, ...]
-
+
def split_gain_il(self, d):
""" Split gain information off AGIPD IL Data
"""
@@ -534,9 +563,9 @@ class AgipdCorrections:
cme, cms = hist_ends(casic[idxm])
asic = d[(i + 1) * 64 - rsize:(i + 1) * 64 - 1,
- 64 - rsize:64].flatten()
+ 64 - rsize:64].flatten()
asic_g = g[(i + 1) * 64 - rsize:(i + 1) * 64 - 1,
- 64 - rsize:64].flatten()
+ 64 - rsize:64].flatten()
idxm = asic_g == 1
me, ms = hist_ends(asic[idxm])
@@ -648,9 +677,9 @@ class AgipdCorrections:
for j in range(2):
fidx = gain[k, i * 64:(i + 1) * 64,
- j * 64:(j + 1) * 64] == 1
+ j * 64:(j + 1) * 64] == 1
fidx_h = gain[k, i * 64:(i + 1) * 64,
- j * 64:(j + 1) * 64] == 0
+ j * 64:(j + 1) * 64] == 0
# need at least two pixels in medium gain to be able to
# cluster in two groups
@@ -705,7 +734,7 @@ class AgipdCorrections:
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]
+ j * 64:(j + 1) * 64]
mim = asic[fidx]
asicb = bg_h
mimb = asicb[fidx]
@@ -747,24 +776,29 @@ class AgipdCorrections:
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
+ j * 64:(j + 1) * 64] = asic_msk
return im, gain, snow_mask
-
+
def make_noisy_adc_mask(self, bmask):
""" Mask entire ADC if they are noise above a relative threshold
- """
- msk = np.count_nonzero(((bmask & BadPixels.NOISE_OUT_OF_THRESHOLD.value != 0)
- | (bmask & BadPixels.OFFSET_NOISE_EVAL_ERROR.value != 0)).astype(np.int) , axis=0)
+ """
+ msk = np.count_nonzero(
+ ((bmask & BadPixels.NOISE_OUT_OF_THRESHOLD.value != 0)
+ | (bmask & BadPixels.OFFSET_NOISE_EVAL_ERROR.value != 0)).astype(
+ np.int), axis=0)
fmask = np.zeros_like(msk).astype(np.uint32)
- adc_i = bmask.shape[1]//8
- adc_j = bmask.shape[2]//8
+ adc_i = bmask.shape[1] // 8
+ adc_j = bmask.shape[2] // 8
for i in range(adc_i):
for j in range(adc_j):
- adc_cnt = np.count_nonzero(msk[i*adc_i:(i+1)*adc_i, j*adc_j:(j+1)*adc_j] > bmask.shape[0]//32)
- if adc_cnt / (adc_i*adc_j) >= self.mask_noisy_adc:
- fmask[i*adc_i:(i+1)*adc_i, j*adc_j:(j+1)*adc_j] = BadPixels.NOISY_ADC.value
- return fmask
+ adc_cnt = np.count_nonzero(
+ msk[i * adc_i:(i + 1) * adc_i, j * adc_j:(j + 1) * adc_j] >
+ bmask.shape[0] // 32)
+ if adc_cnt / (adc_i * adc_j) >= self.mask_noisy_adc:
+ fmask[i * adc_i:(i + 1) * adc_i,
+ j * adc_j:(j + 1) * adc_j] = BadPixels.NOISY_ADC.value
+ return fmask
def correct_agipd(self, irange):
""" Correct Raw AGIPD data for offset and relative gain effects
@@ -788,10 +822,12 @@ class AgipdCorrections:
status = np.squeeze(
self.infile[agipd_base + "/image/status"][irange, ...])
cellid = np.squeeze(
- np.array(self.infile[agipd_base + "/image/cellId"][irange, ...]))
+ np.array(
+ self.infile[agipd_base + "/image/cellId"][irange, ...]))
length = np.squeeze(
- np.array(self.infile[agipd_base + "/image/length"][irange, ...]))
-
+ np.array(
+ self.infile[agipd_base + "/image/length"][irange, ...]))
+
# this far end of the image index range we are working on
nidx = int(cidx + irange.size)
else:
@@ -819,7 +855,7 @@ class AgipdCorrections:
status = status[0::2]
cellid = cellid[0::2]
length = length[0::2]
-
+
# first evaluate the gain into 0, 1, 2 --> high, medium, low
gain = np.zeros(ga.shape, np.uint8)
gain[...] = 0
@@ -828,18 +864,20 @@ class AgipdCorrections:
# 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
-
+ gain[(ga > t1[cellid, ...]) & (
+ t1[cellid, ...] > 1.05 * t0[cellid, ...])] = 2
+
# force into high or medium gain if requested
if self.force_mg_if_below is not None and self.force_mg_if_below > 0:
- gain[(gain == 2) & ((im - self.offset[cellid, ...,1]) < self.force_mg_if_below)] = 1
-
+ gain[(gain == 2) & ((im - self.offset[
+ cellid, ..., 1]) < self.force_mg_if_below)] = 1
+
if self.force_hg_if_below is not None and self.force_hg_if_below > 0:
- gain[(im - self.offset[cellid, ...,0]) < self.force_hg_if_below] = 0
-
+ gain[(gain > 0) & ((im - self.offset[
+ cellid, ..., 0]) < self.force_hg_if_below)] = 0
+
for i in range(3):
self.gain_stats[i] = np.sum(gain == i)
-
# check if any data has zero standard deviation, and mark this in
# the bad pixel maks
@@ -930,6 +968,10 @@ class AgipdCorrections:
if self.do_rel_gain:
im *= rel_cor
+ if self.adjust_mg_baseline:
+ mgbc = self.mg_bl_adjust[cellid, ...]
+ im[gain == 1] += 1.5 * mgbc[gain == 1]
+
# try to identify snowy pixels at this point
if self.melt_snow is not False:
ms = self.melt_snow
@@ -987,14 +1029,14 @@ class AgipdCorrections:
if self.adc_mask is None:
self.adc_mask = self.make_noisy_adc_mask(msk)
msk |= self.adc_mask
-
+
# now write out the data
sd = 1 if not self.il_mode else 2
-
+
if not self.karabo_data_mode:
- self.ddset[cidx//sd:nidx//sd, ...] = im
- self.gdset[cidx//sd:nidx//sd, ...] = gain
- self.mdset[cidx//sd:nidx//sd, ...] = msk
+ self.ddset[cidx // sd:nidx // sd, ...] = im
+ self.gdset[cidx // sd:nidx // sd, ...] = gain
+ self.mdset[cidx // sd:nidx // sd, ...] = msk
self.outfile[agipd_base + "image/cellId"][cidx:nidx] = cellid
self.outfile[agipd_base + "image/trainId"][cidx:nidx] = trainId
@@ -1043,7 +1085,7 @@ class AgipdCorrections:
valid = status != 0
last_index = int(last[status != 0][-1]) + 1
first_index = int(first[status != 0][0])
-
+
idxtrains = np.squeeze(self.infile["/INDEX/trainId"])
medianTrain = np.nanmedian(idxtrains)
lowok = (idxtrains > medianTrain - 1e4)
@@ -1068,7 +1110,7 @@ class AgipdCorrections:
agipd_base = self.agipd_base
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)
@@ -1078,8 +1120,8 @@ class AgipdCorrections:
allcells = allcells[can_calibrate]
firange = np.arange(first_index, last_index)
firange = firange[can_calibrate]
-
- self.oshape = (firange.size if not self.il_mode else firange.size//2,
+
+ self.oshape = (firange.size if not self.il_mode else firange.size // 2,
single_image.shape[1],
single_image.shape[2])
self.firange = firange
@@ -1099,7 +1141,8 @@ class AgipdCorrections:
if self.il_mode:
firange = firange[0::2]
alltrains = self.infile[agipd_base + "image/trainId"]
- alltrains = np.squeeze(alltrains[first_index:last_index, ...])
+ 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",
@@ -1132,15 +1175,33 @@ class AgipdCorrections:
# sanitize indices
for do in ["image", ]:
- print(alltrains.shape, firange.max())
- uq, fidxv, cntsv = np.unique(alltrains[firange-firange[0]], return_index=True,
- return_counts=True)
- fidx = np.zeros(self.valid.shape, fidxv.dtype)
- fidx[self.valid] = fidxv
-
- cnts = np.zeros(self.valid.shape, cntsv.dtype)
- cnts[self.valid] = cntsv
-
+
+ uq, fidxv, cntsv = np.unique(alltrains[firange - firange[0]],
+ return_index=True,
+ return_counts=True)
+
+ duq = (uq[1:] - uq[:-1]).astype(np.int64)
+
+ cfidxv = [fidxv[0], ]
+ ccntsv = [cntsv[0], ]
+ for i, du in enumerate(duq.tolist()):
+ if du > 1000:
+ du = 1
+ cntsv[i] = 0
+ cfidxv += [0] * (du - 1) + [fidxv[i + 1], ]
+ ccntsv += [0] * (du - 1) + [cntsv[i + 1], ]
+
+ mv = len(cfidxv)
+ fidx = np.zeros(len(cfidxv), fidxv.dtype)
+ fidx[self.valid[:mv]] = np.array(cfidxv)[self.valid[:mv]]
+
+ for i in range(len(fidx) - 1, 2, -1):
+ if fidx[i - 1] == 0 and fidx[i] != 0:
+ fidx[i - 1] = fidx[i]
+
+ cnts = np.zeros(len(cfidxv), cntsv.dtype)
+ cnts[self.valid[:mv]] = np.array(ccntsv)[self.valid[:mv]]
+
self.outfile.create_dataset(idx_base + "{}/first".format(do),
fidx.shape,
dtype=fidx.dtype,
@@ -1181,8 +1242,8 @@ class AgipdCorrections:
fsz = list(firange.shape)
if self.il_mode:
- fsz[0] = fsz[0]//2
-
+ fsz[0] = fsz[0] // 2
+
self.outfile.create_dataset(agipd_base + "image/cellId", fsz,
dtype=np.uint16, fletcher32=True)
self.outfile.create_dataset(agipd_base + "image/trainId", fsz,
@@ -1311,7 +1372,7 @@ class AgipdCorrections:
(cal_db_interface, creation_time, max_cells_db,
bias_voltage, photon_energy) = dbparms
max_cells_db_dark = max_cells_db
- timeout_cal_db = 300000
+ timeout_cal_db = 30000
elif len(dbparms) == 6:
(cal_db_interface, creation_time, max_cells_db,
@@ -1319,7 +1380,7 @@ class AgipdCorrections:
if max_cells_db_dark is None:
max_cells_db_dark = max_cells_db
- timeout_cal_db = 300000
+ timeout_cal_db = 30000
else:
(cal_db_interface, creation_time, max_cells_db,
bias_voltage, photon_energy, max_cells_db_dark,
@@ -1328,22 +1389,24 @@ class AgipdCorrections:
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)))])
+ # 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)))])
dinstance = getattr(Detectors, self.cal_det_instance)
- offset = get_constant_from_db(getattr(dinstance, 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)),
- cal_db_interface,
- creation_time=creation_time,
- timeout=timeout_cal_db)
+ offset, when = get_constant_from_db_and_time(getattr(dinstance, 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)),
+ cal_db_interface,
+ creation_time=creation_time,
+ timeout=timeout_cal_db)
noise = get_constant_from_db(getattr(dinstance, qm),
Constants.AGIPD.Noise(),
@@ -1379,8 +1442,9 @@ class AgipdCorrections:
# 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
+ self.initialize(offset, None, bpixels, noise, thresholds, None,
+ slopesPC=None)
+ return when
# load also non-dark constants from the database
slopesPC = get_constant_from_db(getattr(dinstance, qm),
@@ -1408,6 +1472,7 @@ class AgipdCorrections:
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,
@@ -1436,9 +1501,11 @@ class AgipdCorrections:
if slopesFF.shape[2] > 32:
xraygain = np.nanmedian(slopesFF[..., 32:min(96, self.max_cells)],
axis=2)
- else:
+ elif slopesFF.shape[2] > 2:
xraygain = np.nanmedian(slopesFF[..., :min(96, self.max_cells)],
axis=2)
+ else:
+ xraygain = np.squeeze(slopesFF[..., 0])
# relative X-ray gain is then each pixels gain, by the median gain
# of all pixels
@@ -1486,7 +1553,8 @@ class AgipdCorrections:
photon_energy,
pixels_x=512,
pixels_y=128,
- beam_energy=None), # noqa
+ beam_energy=None),
+ # noqa
np.zeros((128, 512, max_cells_db)),
cal_db_interface,
creation_time=creation_time,
@@ -1495,7 +1563,9 @@ class AgipdCorrections:
bpixels |= bpff.astype(np.uint32)[..., :bpixels.shape[2], None]
self.initialize(offset, rel_gain, bpixels, noise, thresholds,
- base_offset)
+ base_offset, slopesPC=slopesPC)
+
+ return when
def initialize_from_file(self, filename, qm, with_dark=True):
""" Initialize calibration constants from a calibration file
@@ -1570,13 +1640,17 @@ class AgipdCorrections:
with h5py.File(filename, "r") as calfile:
bpixels = calfile["{}/{}/data".format(qm, "BadPixelsFF")][()]
- bpixels |= np.moveaxis(
- calfile["{}/{}/data".format(qm, "BadPixelsPC")][()], 0, 2)
+ try:
+ bpixels |= np.moveaxis(
+ calfile["{}/{}/data".format(qm, "BadPixelsPC")][()], 0, 2)
+ except:
+ pass
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")][()] # noqa
+ thresholds = \
+ calfile["{}/{}/data".format(qm, "ThresholdsDark")][()] # noqa
bpixels |= calfile["{}/{}/data".format(qm, "BadPixelsDark")][
()]
@@ -1597,9 +1671,11 @@ class AgipdCorrections:
slopesFF = slopesFF[..., 0]
if slopesFF.shape[2] > 32:
xraygain = np.nanmedian(slopesFF[..., 32:], axis=2)
- else:
+ elif slopesFF.shape[2] > 2:
xraygain = np.nanmedian(
slopesFF[..., :min(96, self.max_cells)], axis=2)
+ else:
+ xraygain = np.squeeze(slopesFF[..., 0])
xraygain /= np.nanmedian(xraygain)
pcrel = pc_high_m / np.nanmedian(pc_high_m, axis=2)[..., None]
@@ -1614,4 +1690,4 @@ class AgipdCorrections:
slopesPC[..., :self.max_cells, 4], xrayOffset]
self.initialize(offsets, rel_gain, bpixels, noises, thresholds,
- base_offset)
+ base_offset, slopesPC=slopesPC)
diff --git a/cal_tools/cal_tools/lpdlib.py b/cal_tools/cal_tools/lpdlib.py
index 5c1a0fb6b8f390400c3fa9993227aca3cc72128b..d0c557aab8964aa847b3aba13d2b063006062254 100644
--- a/cal_tools/cal_tools/lpdlib.py
+++ b/cal_tools/cal_tools/lpdlib.py
@@ -2,8 +2,8 @@ import copy
import h5py
import numpy as np
-from cal_tools.tools import get_constant_from_db
from cal_tools.enums import BadPixels
+from cal_tools.tools import get_constant_from_db, get_constant_from_db_and_time
from iCalibrationDB import Constants, Conditions, Detectors
@@ -41,7 +41,7 @@ class LpdCorrections:
raw_fmt_version=2, chunk_size=512,
h5_data_path="INSTRUMENT/FXE_DET_LPD1M-1/DET/{}CH0:xtdf/",
h5_index_path="INDEX/FXE_DET_LPD1M-1/DET/{}CH0:xtdf/",
- do_ff=True, correct_non_linear=False, karabo_data_mode=False):
+ do_ff=True, correct_non_linear=True, karabo_data_mode=False):
"""
Initialize an LpdCorrections Class
@@ -85,24 +85,19 @@ class LpdCorrections:
self.bins_signal_high_range = bins_signal_high_range
self.cidx = 0
self.do_ff = do_ff
- filter_modules = [11]
+ filter_modules = []
self.filter_cells = [0, 1] if channel in filter_modules else []
- self.cnl = correct_non_linear
+ self.cnl = True # correct_non_linear
self.karabo_data_mode = karabo_data_mode
# emprically determined from APD datasets p900038, r155,r156
- self.cnl_const = {'med':
- {'lin': (0.10960777594700886,
- 859.4778021212404),
- 'repo': (0.11046365063379418,
- 751.5109159229758,
- 1436280.868413923,
- 9966.92701736398)},
- 'high': {'sigm': (1333.6836067255365,
- 8370.936473950766,
- 0.003912416547243528,
- 0.00013306099710934553),
- 'lin': (0.16284255864598732,
- -29.830797125784763)}}
+ # emprically determined from APD datasets p900038, r155,r156
+ self.cnl_const = {
+ 'high': {'A': -0.000815934, 'lam': 0.00811281, 'c': 1908.89,
+ 'm': 0, 'b': 0},
+ 'med': {'A': 0.0999894, 'lam': -0.00137652, 'c': 3107.83,
+ 'm': 3.89982e-06, 'b': -0.116811},
+ 'low': {'A': 0.0119132, 'lam': -0.0002, 'c': 36738.6,
+ 'm': 2.00273e-08, 'b': 0.245537}}
def get_iteration_range(self):
"""Returns a range expression over which to iterate in chunks
@@ -157,7 +152,7 @@ class LpdCorrections:
if not hasattr(self, "mask"):
self.mask = mask
else:
- self.mask |= mask
+ self.mask |= mask[:self.mask.shape[0], ...]
if noise is not None:
self.noise = noise
if flatfield is not None:
@@ -222,7 +217,6 @@ class LpdCorrections:
cells = np.squeeze(irange['image.cellId'])
length = np.squeeze(irange['image.length'])
-
# split gain and image info into separate arrays
im, gain = self.split_gain(im[:, 0, ...])
@@ -269,35 +263,30 @@ class LpdCorrections:
im /= self.flatfield[None, :, :]
# hacky way of smoothening transition region between med and low
- #im[gain == 2] -= im[gain == 2] * cfac[gain == 2]
-
+ # im[gain == 2] -= im[gain == 2] * cfac[gain == 2]
+
# perform non-linear corrections if requested
if self.cnl:
- largs = self.cnl_const['high']['lin']
- sargs = self.cnl_const['high']['sigm']
- n, c, b, k = sargs
- y = im[(gain==0) & (im > 1000)]
- x = (-np.arctanh(b-y/n+1)+c*k)/k
- dy = largs[0]*x+largs[1] - n*(1+b+np.tanh((x-c)*k))
- im[(gain==0) & (im > 1000)] += dy
-
- mhg = largs[0]
- bhg = largs[1]
-
- largs = self.cnl_const['med']['lin']
- rargs = self.cnl_const['med']['repo']
- c, k, j, n = rargs
- j *= -1
- mf = mhg/largs[0]
-
- y = im[(gain==1)]
- x = (np.sqrt(c**2*n**2+4*c*j+2*c*n*(k-y)+(k-y)**2)+c*n-k+y)/(2*c) - (bhg-largs[1])/largs[0]
- dy = (mf*largs[0]-c)*x + largs[1] - k + j/(x-n)
-
- im[(gain==1)] += (bhg-largs[1]) + dy
-
- #im[(gain==1) & (im >= 3500)] *= mf
- #im[(gain==1) & (im >= 3500)] += (bhg-largs[1])
+ def lin_exp_fun(x, m, b, A, lam, c):
+ return m * x + b + A * np.exp(lam * (x - c))
+
+ x = im[(gain == 0)]
+ cnl = self.cnl_const['high']
+ cf = lin_exp_fun(x, cnl['m'], cnl['b'], cnl['A'], cnl['lam'],
+ cnl['c'])
+ im[(gain == 0)] -= np.maximum(cf, -0.2) * x
+
+ x = im[(gain == 1)]
+ cnl = self.cnl_const['med']
+ cf = lin_exp_fun(x, cnl['m'], cnl['b'], cnl['A'], cnl['lam'],
+ cnl['c'])
+ im[(gain == 1)] -= np.minimum(cf, 0.2) * x
+
+ x = im[(gain == 2)]
+ cnl = self.cnl_const['low']
+ cf = lin_exp_fun(x, cnl['m'], cnl['b'], cnl['A'], cnl['lam'],
+ cnl['c'])
+ im[(gain == 2)] -= np.minimum(cf, 0.45) * x
# create bad pixels masks, here non-finite values
bidx = ~np.isfinite(im)
@@ -378,11 +367,18 @@ class LpdCorrections:
status = np.squeeze(self.infile[lpd_base + "image/status"])
if np.count_nonzero(status != 0) == 0:
raise IOError("File {} has no valid counts".format(
- self.infile))
+ self.infile))
last = np.squeeze(self.infile[lpd_base + "image/last"])
valid = status != 0
- last_index = int(last[status != 0][-1])
- first_index = int(last[status != 0][0])
+
+ idxtrains = np.squeeze(self.infile["/INDEX/trainId"])
+ medianTrain = np.nanmedian(idxtrains)
+ lowok = (idxtrains > medianTrain - 1e4)
+ highok = (idxtrains < medianTrain + 1e4)
+ valid &= lowok & highok
+
+ last_index = int(last[valid][-1])
+ first_index = int(last[valid][0])
else:
raise AttributeError(
"Not a known raw format version: {}".format(self.index_v))
@@ -433,11 +429,6 @@ class LpdCorrections:
dont_copy = [lpd_base + "image/{}".format(do)
for do in dont_copy]
- # dont_copy += [idx_base+"{}/first".format(do)
- # for do in ["detector", "header", "image", "trailer"]]
- # dont_copy += [idx_base+"{}/count".format(do)
- # for do in ["detector", "header", "image", "trailer"]]
-
dont_copy += [idx_base + "{}/first".format(do)
for do in ["image"]]
dont_copy += [idx_base + "{}/count".format(do)
@@ -457,8 +448,32 @@ class LpdCorrections:
self.infile.visititems(visitor)
# sanitize indices
for do in ["image", ]:
- uq, fidx, cnts = np.unique(alltrains[firange], return_index=True,
- return_counts=True)
+ uq, fidxv, cntsv = np.unique(alltrains[firange - firange[0]],
+ return_index=True,
+ return_counts=True)
+
+ duq = (uq[1:] - uq[:-1]).astype(np.int64)
+
+ cfidxv = [fidxv[0], ]
+ ccntsv = [cntsv[0], ]
+ for i, du in enumerate(duq.tolist()):
+ if du > 1000:
+ du = 1
+ cntsv[i] = 0
+ cfidxv += [0] * (du - 1) + [fidxv[i + 1], ]
+ ccntsv += [0] * (du - 1) + [cntsv[i + 1], ]
+
+ mv = len(cfidxv)
+ fidx = np.zeros(len(cfidxv), fidxv.dtype)
+ fidx[self.valid[:mv]] = np.array(cfidxv)[self.valid[:mv]]
+
+ for i in range(len(fidx) - 1, 2, -1):
+ if fidx[i - 1] == 0 and fidx[i] != 0:
+ fidx[i - 1] = fidx[i]
+
+ cnts = np.zeros(len(cfidxv), cntsv.dtype)
+ cnts[self.valid[:mv]] = np.array(ccntsv)[self.valid[:mv]]
+
self.outfile.create_dataset(idx_base + "{}/first".format(do),
fidx.shape,
dtype=fidx.dtype,
@@ -585,23 +600,17 @@ class LpdCorrections:
(cal_db_interface, creation_time, max_cells_db, capacitor,
bias_voltage, photon_energy, timeout) = dbparms
- 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)))])
- print("Connecting to interface: {}".format(cal_db_interface))
-
- offset = get_constant_from_db(getattr(Detectors.LPD1M1, qm),
- Constants.LPD.Offset(),
- Conditions.Dark.LPD(
- memory_cells=max_cells_db,
- bias_voltage=bias_voltage,
- capacitor=capacitor),
- np.zeros((256, 256, max_cells_db, 3)),
- cal_db_interface,
- creation_time=creation_time,
- timeout=timeout)
+ offset, when = get_constant_from_db_and_time(
+ getattr(Detectors.LPD1M1, qm),
+ Constants.LPD.Offset(),
+ Conditions.Dark.LPD(
+ memory_cells=max_cells_db,
+ bias_voltage=bias_voltage,
+ capacitor=capacitor),
+ np.zeros((256, 256, max_cells_db, 3)),
+ cal_db_interface,
+ creation_time=creation_time,
+ timeout=timeout)
noise = get_constant_from_db(getattr(Detectors.LPD1M1, qm),
Constants.LPD.Noise(),
@@ -628,7 +637,7 @@ class LpdCorrections:
# done loading constants and returning
if only_dark:
self.initialize(offset, None, None, bpixels, noise, None)
- return
+ return when
slopesCI = get_constant_from_db(getattr(Detectors.LPD1M1, qm),
Constants.LPD.SlopesCI(),
@@ -685,6 +694,7 @@ class LpdCorrections:
bpixels |= bpix[..., None]
self.initialize(offset, rel_gains, rel_gain_offset, bpixels, noise,
flat_fields)
+ return when
def initialize_from_file(self, filename, qm, with_dark=True):
""" Initialize calibration constants from a calibration file
diff --git a/cal_tools/cal_tools/tools.py b/cal_tools/cal_tools/tools.py
index 2e8961c002d7685df06ef97a0671f0ad6db4ece2..e1e68c1589b565d7586f2e575b6ebd5e92d69cb3 100644
--- a/cal_tools/cal_tools/tools.py
+++ b/cal_tools/cal_tools/tools.py
@@ -9,9 +9,8 @@ from subprocess import CalledProcessError, check_call, check_output
from textwrap import dedent
from time import sleep
-from jinja2 import Template
import tabulate
-
+from jinja2 import Template
from xfel_calibrate import settings
@@ -38,24 +37,27 @@ def combine_report(run_path, calibration):
makedirs(sphinx_path)
direntries = listdir(run_path)
direntries.sort(key=natural_keys)
-
+
for entry in direntries:
-
+
if isfile("{}/{}".format(run_path, entry)):
name, ext = splitext("{}".format(entry))
-
+
if ext == ".rst":
- comps = name.split("__")
+ comps = name.split("__")
if len(comps) >= 3:
- group, name_param, conc_param = "_".join(comps[:-2]), comps[-2], comps[-1]
+ group, name_param, conc_param = "_".join(comps[:-2]), \
+ comps[-2], comps[-1]
else:
- group, name_param, conc_param = comps[0], "None", "None"
+ group, name_param, conc_param = comps[0], "None", "None"
- with open("{}/{}.rst".format(sphinx_path, group), "a") as gfile:
+ with open("{}/{}.rst".format(sphinx_path, group),
+ "a") as gfile:
if conc_param != "None":
- title = "{}. {} = {}".format(calibration, name_param, conc_param)
+ title = "{}. {} = {}".format(calibration, name_param,
+ conc_param)
gfile.write(title + "\n")
- gfile.write( "=" *len (title) + "\n")
+ gfile.write("=" * len(title) + "\n")
gfile.write("\n")
with open("{}/{}".format(run_path, entry), "r") as ifile:
skip_next = False
@@ -65,14 +67,16 @@ def combine_report(run_path, calibration):
continue
if conc_param == "None":
gfile.write(line)
- elif " ".join(calibration.split()) != " ".join(line.split()):
+ elif " ".join(calibration.split()) != " ".join(
+ line.split()):
gfile.write(line)
else:
skip_next = True
-
+
gfile.write("\n\n")
if isdir("{}/{}".format(run_path, entry)):
- copytree("{}/{}".format(run_path, entry), "{}/{}".format(sphinx_path, entry))
+ copytree("{}/{}".format(run_path, entry),
+ "{}/{}".format(sphinx_path, entry))
return sphinx_path
@@ -268,7 +272,7 @@ def make_report(run_path, tmp_path, out_path, project, author, version,
print("Removing temporary files at: {}".format(tmp_path))
rmtree(tmp_path)
-
+
def make_titlepage(sphinx_path, project, data_path):
"""
Create title page for report using template
@@ -288,7 +292,6 @@ def make_titlepage(sphinx_path, project, data_path):
def finalize(joblist, finaljob, run_path, out_path, project, calibration,
author, version, report_to, data_path='Unknown'):
-
print("Waiting on jobs to finish: {}".format(joblist))
while True:
found_jobs = set()
@@ -302,7 +305,7 @@ def finalize(joblist, finaljob, run_path, out_path, project, calibration,
sleep(10)
prepare_plots(run_path)
- make_timing_summary(run_path, joblist+[str(finaljob)])
+ make_timing_summary(run_path, joblist + [str(finaljob)])
sphinx_path = combine_report(run_path, calibration)
make_titlepage(sphinx_path, project, data_path)
make_report(sphinx_path, run_path, out_path, project, author, version,
@@ -321,13 +324,13 @@ def parse_runs(runs, return_type=str):
elif isinstance(runs, (list, tuple)):
pruns = runs
else:
- pruns = [runs,]
-
+ pruns = [runs, ]
+
if return_type is str:
return ["r{:04d}".format(r) for r in pruns]
else:
return pruns
-
+
def run_prop_seq_from_path(filename):
run = re.findall(r'.*r([0-9]{4}).*', filename)
@@ -356,11 +359,12 @@ def map_modules_from_files(filelist, file_inset, quadrants, modules_per_quad):
module_files[name].put(file)
total_sequences += 1
total_file_size += path.getsize(file)
-
+
return module_files, mod_ids, total_sequences, total_file_size
-def gain_map_files(in_folder, runs, sequences, file_inset, quadrants, modules_per_quad):
+def gain_map_files(in_folder, runs, sequences, file_inset, quadrants,
+ modules_per_quad):
total_sequences = 0
total_file_size = 0
gain_mapped_files = OrderedDict()
@@ -369,7 +373,7 @@ def gain_map_files(in_folder, runs, sequences, file_inset, quadrants, modules_pe
dirlist = listdir(ginfolder)
file_list = []
for entry in dirlist:
- #only h5 file
+ # only h5 file
abs_entry = "{}/{}".format(ginfolder, entry)
if path.isfile(abs_entry) and path.splitext(abs_entry)[1] == ".h5":
@@ -379,13 +383,15 @@ def gain_map_files(in_folder, runs, sequences, file_inset, quadrants, modules_pe
for seq in sequences:
if "{:05d}.h5".format(seq) in abs_entry:
file_list.append(path.abspath(abs_entry))
-
- mapped_files, mod_ids, seq, fs = map_modules_from_files(file_list, file_inset,
- quadrants, modules_per_quad)
+
+ mapped_files, mod_ids, seq, fs = map_modules_from_files(file_list,
+ file_inset,
+ quadrants,
+ modules_per_quad)
total_sequences += seq
total_file_size += fs
gain_mapped_files[gain] = mapped_files
- return gain_mapped_files, total_sequences, total_file_size/1e9
+ return gain_mapped_files, total_sequences, total_file_size / 1e9
# from https://github.com/jupyter/notebook/issues/1000
@@ -405,6 +411,7 @@ try: # Python 3
except ImportError: # Python 2
import warnings
from IPython.utils.shimmodule import ShimWarning
+
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=ShimWarning)
from IPython.html.notebookapp import list_running_servers
@@ -427,28 +434,38 @@ def get_notebook_name():
return relative_path
except:
return os.environ.get("CAL_NOTEBOOK_NAME", "Unknown Notebook")
-
-
+
+
def get_dir_creation_date(directory, run):
import os
import datetime
- creation_time = os.stat("{}/r{:04d}".format(directory, run)).st_ctime
+ creation_time = os.stat("{}/r{:04d}".format(directory, run)).st_ctime
creation_time = datetime.datetime.fromtimestamp(creation_time)
- return(creation_time)
+ return creation_time
def get_dir_creation_date(directory, run):
import os
import datetime
- creation_time = os.stat("{}/r{:04d}".format(directory, run)).st_mtime
+ creation_time = os.stat("{}/r{:04d}".format(directory, run)).st_mtime
creation_time = datetime.datetime.fromtimestamp(creation_time)
- return(creation_time)
+ return creation_time
+
+
+def get_random_db_interface(cal_db_interface):
+ if "#" in cal_db_interface:
+ prot, serv, ran = cal_db_interface.split(":")
+ r1, r2 = ran.split("#")
+ return ":".join([prot, serv, str(np.random.randint(int(r1), int(r2)))])
+ return cal_db_interface
already_printed = {}
+
+
def get_from_db(device, constant, condition, empty_constant,
- cal_db_interface, creation_time = None,
- verbosity=1, timeout=30000):
+ cal_db_interface, creation_time=None,
+ verbosity=1, timeout=30000, ntries=120):
"""
Return calibration constants and metadata requested from CalDB
@@ -460,32 +477,53 @@ def get_from_db(device, constant, condition, empty_constant,
:param creation_time: Latest time for constant to be created
:param verbosity: Level of verbosity (0 - silent)
:param timeout: Timeout for zmq request
+ :param ntries: number of tries to contact the database
:return: Calibration constant, metadata
"""
- from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
-
+ from iCalibrationDB import ConstantMetaData, Versions
+
if device:
- metadata = ConstantMetaData()
+ metadata = ConstantMetaData()
metadata.calibration_constant = constant
metadata.detector_condition = condition
if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device)
+ while ntries > 0:
+ this_interface = get_random_db_interface(cal_db_interface)
+ try:
+ metadata.retrieve(this_interface, timeout=6000000,
+ meta_only=True)
+ break
+ except Exception as e:
+ if verbosity > 0:
+ print(e)
+ ntries -= 1
else:
- metadata.calibration_constant_version = Versions.Timespan(device=device,
- start=creation_time)
- creation_time = creation_time.isoformat()
-
- try:
- metadata.retrieve(cal_db_interface, when=creation_time, timeout=timeout)
+ metadata.calibration_constant_version = Versions.Timespan(
+ device=device,
+ start=creation_time)
+
+ while ntries > 0:
+ this_interface = get_random_db_interface(cal_db_interface)
+ try:
+ metadata.retrieve(this_interface,
+ when=creation_time.isoformat(),
+ timeout=6000000, meta_only=True)
+ break
+ except Exception as e:
+ if verbosity > 0:
+ print(e)
+
+ ntries -= 1
+ if ntries > 0:
if verbosity > 0:
if constant.name not in already_printed or verbosity > 1:
already_printed[constant.name] = True
+ begin_at = metadata.calibration_constant_version.begin_at
print("{} was injected on: {}".format(constant.name,
- metadata.calibration_constant_version.begin_at))
+ begin_at))
return constant.data, metadata
- except Exception as e:
- if verbosity > 0:
- print(e)
+ else:
return empty_constant, metadata
else:
return empty_constant, None
@@ -493,17 +531,33 @@ def get_from_db(device, constant, condition, empty_constant,
def get_constant_from_db(device, constant, condition, empty_constant,
cal_db_interface, creation_time=None,
- print_once=True, timeout=30000):
+ print_once=True, timeout=30000, ntries=120):
"""
Return calibration constants requested from CalDB
"""
data, _ = get_from_db(device, constant, condition, empty_constant,
- cal_db_interface, creation_time,
- int(print_once), timeout)
+ cal_db_interface, creation_time,
+ int(print_once), timeout, ntries)
return data
+def get_constant_from_db_and_time(device, constant, condition, empty_constant,
+ cal_db_interface, creation_time=None,
+ print_once=True, timeout=30000, ntries=120):
+ """
+ Return calibration constants requested from CalDB, alongside injection
+ time
+ """
+ data, m = get_from_db(device, constant, condition, empty_constant,
+ cal_db_interface, creation_time,
+ int(print_once), timeout, ntries)
+ if m:
+ return data, m.calibration_constant_version.begin_at
+ else:
+ return data, None
+
+
def tex_escape(text):
"""
Escape latex special characters found in the text