Skip to content
Snippets Groups Projects
Commit 4370a8ea authored by Steffen Hauf's avatar Steffen Hauf
Browse files

Resolved regressions

parent b5b55ae7
No related branches found
No related tags found
1 merge request!109Add a manual slow control data to Jungfrau notebooks
Hide whitespace changes
Inline Side-by-side
notebooks/Jungfrau/Jungfrau_Gain_Correct_and_Verify_NBC.ipynb
+ 9
− 10
View file @ 4370a8ea
%% Cell type:markdown id: tags:
# Jungfrau Offline Correction #
Author: European XFEL Detector Group, Version: 0.1
Offline Calibration for the Jungfrau Detector
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/FXE/201802/p002271/ra" # the folder to read data from, required
run = 132 # runs to process, required
# out_folder = "/gpfs/exfel/data/scratch/haufs/test/" # the folder to output to, required
out_folder = "/gpfs/exfel/data/scratch/xcal/test/"
out_folder = "/gpfs/exfel/data/scratch/haufs/test/" # the folder to output to, required
calfile = "" # path to calibration file. Leave empty if all data should come from DB, not actually used, makes automode happy
sequences = [-1] # sequences to correct, set to -1 for all, range allowed
mem_cells = 1 # memory cells in data, not actually used, makes automode happy
overwrite = True # set to True if existing data should be overwritten
no_relative_gain = False # do not do relative gain correction
cluster_profile = "noDB"
cluster_profile = "noDB" # cluster profile to use
bias_voltage = 90 # will be overwritten by value in file
cal_db_interface = "tcp://max-exfl016:8016" #"tcp://max-exfl016:8015#8025" # the database interface to use
use_dir_creation_date = True # use the creation data of the input dir for database queries
sequences_per_node = 5 # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel
photon_energy = 9.2 # photon energy in keV
nodb = False # if set only file-based constants will be used
chunk_size_idim = 1 # chunking size of imaging dimension, adjust if user software is sensitive to this.
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # template to use for file name, double escape sequence number
path_template_control = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # template to use for file name, double escape sequence number
integration_time = 1000 # integration time in us, will be overwritten by value in file
h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data' # path in H5 file under which images are located
gmapfile = "/gpfs/exfel/data/scratch/xcal/jfgain/gainMaps_M233.h5" #temporary gain calibration file, not in the DB; leave empty if using DB
memcells = 1
karabo_id = "FXE_XAD_JF1M"
receiver_id = "RECEIVER"
control_id = "CONTROL"
db_module = "Jungfrau_M233"
path_inset = "JNGFR02"
path_inset_control = "JNGFR01"
memcells = 1 # number of memory cells
karabo_id = "FXE_XAD_JF1M" # karabo prefix of Jungfrau devices
receiver_id = "RECEIVER" # inset for receiver devices
control_id = "CONTROL" # inset for control devices
db_module = "Jungfrau_M233" # ID of module in calibration database
path_inset = "JNGFR02" # file inset for image data
path_inset_control = "JNGFR01" # file inset for control data
manual_slow_data = False # set this flag to not use slow control data from file, but manual entries
def balance_sequences(in_folder, run, sequences, sequences_per_node):
import glob
import re
import numpy as np
if sequences_per_node != 0:
sequence_files = glob.glob("{}/r{:04d}/*-S*.h5".format(in_folder, run))
seq_nums = set()
for sf in sequence_files:
seqnum = re.findall(r".*-S([0-9]*).h5", sf)[0]
seq_nums.add(int(seqnum))
seq_nums -= set(sequences)
return [l.tolist() for l in np.array_split(list(seq_nums),
len(seq_nums)//sequences_per_node+1)]
else:
return sequences
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
import warnings
warnings.filterwarnings('ignore')
h5path = h5path.format(karabo_id, receiver_id)
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import XFELDetAna.xfelpyanatools as xana
import XFELDetAna.xfelpycaltools as xcal
from cal_tools.enums import BadPixels
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from cal_tools.tools import (get_dir_creation_date)
import os
import h5py
from XFELDetAna.util import env
env.iprofile = cluster_profile
from XFELDetAna.detectors.jungfrau import readerPSI as jfreaderPSI
from XFELDetAna.detectors.jungfrau import reader as jfreader
from XFELDetAna.detectors.jungfrau.jf_chunk_reader import JFChunkReader
from XFELDetAna.detectors.jungfrau.util import non_empty_trains
# so constants relevant for the analysis
cpuCores = 1
is_parallel = False
sensorSize = [1024, 512]
blockSize = [1024, 512]
xRange = [0, 0+sensorSize[0]]
yRange = [0, 0+sensorSize[1]]
gains = [0, 1, 2]
ped_dir = "{}/r{:04d}".format(in_folder, run)
if ped_dir[-1] == "/":
ped_dir = ped_dir[:-1]
out_folder = "{}/{}".format(out_folder, os.path.split(ped_dir)[-1])
if not os.path.exists(out_folder):
os.makedirs(out_folder)
elif not overwrite:
raise AttributeError("Output path exists! Exiting")
fp_name = path_template.format(run, path_inset)
fp_path = '{}/{}'.format(ped_dir, fp_name)
fp_name_contr = path_template_control.format(run, path_inset_control)
fp_path_contr = '{}/{}'.format(ped_dir, fp_name_contr)
filep_size = 500
myRange_P = sequences
path = h5path
if sequences[0] == -1:
sequences = None
print("Reading data from {}".format(fp_path))
print("Run is: {}".format(run))
print("HDF5 path: {}".format(h5path))
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run)
print("Using {} as creation time".format(creation_time))
cal_timeout = 600000 #ms
```
%% Cell type:code id: tags:
``` python
import h5py
if not manual_slow_data:
with h5py.File(fp_path_contr.format(0), 'r') as f:
integration_time = int(f['/RUN/{}/DET/{}/exposureTime/value'.format(karabo_id, control_id)][()]*1e6)
bias_voltage = int(np.squeeze(f['/RUN/{}/DET/{}/vHighVoltage/value'.format(karabo_id, control_id)])[0])
print("Integration time is {} us".format(integration_time))
print("Bias voltage is {} V".format(bias_voltage))
```
%% Cell type:code id: tags:
``` python
dirlist = sorted(os.listdir(ped_dir))
file_list = []
total_sequences = 0
fsequences = []
for entry in dirlist:
#only h5 file
abs_entry = "{}/{}".format(ped_dir, entry)
if os.path.isfile(abs_entry) and os.path.splitext(abs_entry)[1] == ".h5":
if sequences is None:
for seq in range(len(dirlist)):
if path_template.format(run, path_inset).format(seq) in abs_entry:
file_list.append(abs_entry)
total_sequences += 1
fsequences.append(seq)
else:
for seq in sequences:
if path_template.format(run, path_inset).format(seq) in abs_entry:
file_list.append(os.path.abspath(abs_entry))
total_sequences += 1
fsequences.append(seq)
sequences = fsequences
```
%% Cell type:code id: tags:
``` python
import copy
from IPython.display import HTML, display, Markdown, Latex
import tabulate
print("Processing a total of {} sequence files".format(total_sequences))
table = []
for k, f in enumerate(file_list):
table.append((k, f))
md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=["#", "file"])))
```
%% Cell type:code id: tags:
``` python
def get_PSI_gainmaps(fname, dset):
with h5py.File(fname, 'r') as f_g:
gain_map = np.array(f_g[dset])
gain_map = np.transpose(gain_map, (2, 1, 0, 3))
gain_map = np.expand_dims(gain_map, axis=0)
return gain_map
```
%% Cell type:code id: tags:
``` python
## offset
metadata = ConstantMetaData()
offset = Constants.jungfrau.Offset()
metadata.calibration_constant = offset
# set the operating condition
condition = Conditions.Dark.jungfrau(memory_cells=1, bias_voltage=bias_voltage,
integration_time=integration_time)
device = getattr(Detectors, db_module)
metadata.detector_condition = condition
# specify the a version for this constant
if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device)
metadata.retrieve(cal_db_interface)
else:
metadata.calibration_constant_version = Versions.Timespan(device=device,
start=creation_time)
metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)
offset_map = metadata.calibration_constant.data
## noise
metadata = ConstantMetaData()
noise = Constants.jungfrau.Noise()
metadata.calibration_constant = noise
# set the operating condition
condition = Conditions.Dark.jungfrau(memory_cells=1, bias_voltage=bias_voltage,
integration_time=integration_time)
metadata.detector_condition = condition
# specify the a version for this constant
if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device)
metadata.retrieve(cal_db_interface)
else:
metadata.calibration_constant_version = Versions.Timespan(device=device,
start=creation_time)
metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)
noise_map = metadata.calibration_constant.data
## bad pixels
metadata = ConstantMetaData()
bpix = Constants.jungfrau.BadPixelsDark()
metadata.calibration_constant = bpix
# set the operating condition
condition = Conditions.Dark.jungfrau(memory_cells=1, bias_voltage=bias_voltage,
integration_time=integration_time)
metadata.detector_condition = condition
# specify the a version for this constant
if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device)
metadata.retrieve(cal_db_interface)
else:
metadata.calibration_constant_version = Versions.Timespan(device=device,
start=creation_time)
metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)
mask = metadata.calibration_constant.data
```
%% Cell type:code id: tags:
``` python
# gain correction
if gmapfile:
gain_map = get_PSI_gainmaps(gmapfile, 'gain_map_g0')
print('Not Using DB gain maps ')
print('maps from: {}'.format(gmapfile))
```
%% Cell type:code id: tags:
``` python
def copy_and_sanitize_non_cal_data(infile, outfile, h5base):
""" Copy and sanitize data in `infile` that is not touched by `correctLPD`
"""
if h5base.startswith("/"):
h5base = h5base[1:]
dont_copy = ["adc",]
dont_copy = [h5base+"/{}".format(do)
for do in dont_copy]
def visitor(k, item):
if k not in dont_copy:
if isinstance(item, h5py.Group):
outfile.create_group(k)
elif isinstance(item, h5py.Dataset):
group = str(k).split("/")
group = "/".join(group[:-1])
infile.copy(k, outfile[group])
infile.visititems(visitor)
```
%% Cell type:code id: tags:
``` python
import copy
offset_map = np.squeeze(offset_map)
mask = np.squeeze(mask)
fim_data = None
gim_data = None
rim_data = None
msk_data = None
for k, f in enumerate(file_list):
with h5py.File(f, 'r') as infile:
out_file = "{}/{}".format(out_folder, f.split("/")[-1])
out_file = out_file.replace("RAW", "CORR")
with h5py.File(out_file, "w") as ofile:
copy_and_sanitize_non_cal_data(infile, ofile, h5path)
d = infile[h5path+"/adc"][()].astype(np.float32)
if k == 0:
rim_data = np.squeeze(copy.copy(d))
g = infile[h5path+"/gain"][()]
oshape = d.shape
ddset = ofile.create_dataset(h5path+"/adc",
oshape,
chunks=(chunk_size_idim, 1, oshape[1], oshape[2]),
dtype=np.float32)
mskset = ofile.create_dataset(h5path+"/mask",
oshape,
chunks=(chunk_size_idim, 1, oshape[1], oshape[2]),
dtype=np.uint32,
compression="gzip", compression_opts=1, shuffle=True)
g[g==3] = 2
#offset correction
offset = np.choose(g, (offset_map[...,0], offset_map[...,1], offset_map[...,2]))
d -= offset
#gain correction
cal = np.choose(g, (gain_map[..., 0], gain_map[..., 1], gain_map[..., 2]))
d /= cal
ddset[...] = d
msk = np.choose(g, (mask[...,0], mask[...,1], mask[...,2]))
mskset[...] = msk
if k == 0:
fim_data = np.squeeze(copy.copy(d))
gim_data = np.squeeze(copy.copy(g))
msk_data = np.squeeze(copy.copy(msk))
```
%% Cell type:code id: tags:
``` python
def do_2d_plot(data, edges, y_axis, x_axis):
from matplotlib.colors import LogNorm
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
extent = [np.min(edges[1]), np.max(edges[1]),np.min(edges[0]), np.max(edges[0])]
im = ax.imshow(data[::-1,:], extent=extent, aspect="auto", norm=LogNorm(vmin=1, vmax=np.max(data)))
ax.set_xlabel(x_axis)
ax.set_ylabel(y_axis)
cb = fig.colorbar(im)
cb.set_label("Counts")
```
%% Cell type:code id: tags:
``` python
h, ex, ey = np.histogram2d(rim_data.flatten(), gim_data.flatten(),
bins=[100, 4], range=[[0, 10000], [0,4]])
do_2d_plot(h, (ex, ey), "Signal (ADU)", "Gain Bit Value")
```
%% Cell type:markdown id: tags:
### Mean RAW Preview ###
The per pixel mean of the sequence file of RAW data
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.mean(rim_data,axis=0),
vmin=min(0.75*np.median(rim_data[rim_data > 0]), 2000),
vmax=max(1.5*np.median(rim_data[rim_data > 0]), 16000), cmap="jet")
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:markdown id: tags:
### Mean CORRECTED Preview ###
The per pixel mean of the sequence file of CORR data
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.mean(fim_data,axis=0),
vmin=min(0.75*np.median(fim_data[fim_data > 0]), -0.5),
vmax=max(2.*np.median(fim_data[fim_data > 0]), 100), cmap="jet")
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:markdown id: tags:
### Single Train Preview ###
A single image from the first train
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(fim_data[0,...],
vmin=min(0.75*np.median(fim_data[0,...]), -0.5),
vmax=max(2.*np.median(fim_data[0,...]), 1000), cmap="jet")
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:markdown id: tags:
## Signal Distribution ##
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(211)
h = ax.hist(fim_data.flatten(), bins=1000, range=(-100, 1000), log=True)
l = ax.set_xlabel("Signal (keV)")
l = ax.set_ylabel("Counts")
ax = fig.add_subplot(212)
h = ax.hist(fim_data.flatten(), bins=1000, range=(-100, 10000), log=True)
l = ax.set_xlabel("Signal (keV)")
l = ax.set_ylabel("Counts")
```
%% Cell type:markdown id: tags:
### Maximum GAIN Preview ###
The per pixel maximum of the first train of the GAIN data
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.max(gim_data, axis=0), vmin=0,
vmax=3, cmap="jet")
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:markdown id: tags:
## Bad Pixels ##
The mask contains dedicated entries for all pixels and memory cells as well as all three gains stages. Each mask entry is encoded in 32 bits as:
%% Cell type:code id: tags:
``` python
from cal_tools.enums import BadPixels
from IPython.display import HTML, display, Markdown, Latex
import tabulate
table = []
for item in BadPixels:
table.append((item.name, "{:016b}".format(item.value)))
md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=["Bad pixel type", "Bit mask"])))
```
%% Cell type:markdown id: tags:
### Single Image Bad Pixels ###
A single image bad pixel map for the first image of the first train
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
im = ax.imshow(np.log2(msk_data[0,...]), vmin=0, vmax=32, cmap="jet")
cb = fig.colorbar(im, ax=ax)
```
%% Cell type:code id: tags:
``` python
```
......
notebooks/Jungfrau/Jungfrau_dark_analysis_all_gains_NBC.ipynb
+ 7
− 7
View file @ 4370a8ea
%% Cell type:markdown id: tags:
# Jungfrau Dark Image Characterization #
Version: 0.1, Author: M. Ramilli, S. Hauf
Analyzes Jungfrau dark image data to deduce offset, noise and resulting bad pixel maps
%% Cell type:code id: tags:
``` python
in_folder = '/gpfs/exfel/exp/SPB/201921/p002429/raw/' # folder under which runs are located, required
out_folder = '' # path to place reports at, required
sequences = 1 # number of sequence files in that run
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # template to use for file name, double escape sequence number
path_inset = "DA06"
path_inset_control = "DA06"
path_inset = "DA06" # file inset for image data
path_inset_control = "DA06" # file inset for control data
cluster_profile = 'noDB' # the ipcluster profile name
cal_db_interface = 'tcp://max-exfl016:8016' # calibrate db interface to connect to
integration_time = 1000 # integration time in us, will be overwritten by value in file
bias_voltage = 90 # sensor bias voltage in V, will be overwritten by value in file
badpixel_threshold_sigma = 20. # bad pixels defined by values outside n times this std from median
offset_abs_threshold_low = [1000, 10000, 10000] # absolute bad pixel threshold in terms of offset, lower values
offset_abs_threshold_high = [8000, 15000, 15000] # absolute bad pixel threshold in terms of offset, upper values
chunkSize = 10 # iteration chunk size, needs to match or be less than number of images in a sequence file
imageRange = [0, 500] # image range in which to evaluate
memoryCells = 1 # number of memory cells
h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data' # path in H5 file under which images are located
run_high = 0 # run number for G0 dark run, required
run_med = 0 # run number for G1 dark run, required
run_low = 0 # run number for G2 dark run, required
karabo_id = "FXE_XAD_JF500K"
receiver_id = "RECEIVER"
control_id = "CONTROL"
db_module = "Jungfrau_M233"
use_dir_creation_date = True
karabo_id = "FXE_XAD_JF500K" # karabo prefix of Jungfrau devices
receiver_id = "RECEIVER" # inset for receiver devices\
control_id = "CONTROL" # inset for control devices
db_module = "Jungfrau_M233" # ID of module in calibration database
use_dir_creation_date = True # use dir creation date
manual_slow_data = False # set this flag to not use slow control data from file, but manual entries
```
%% Cell type:code id: tags:
``` python
import warnings
warnings.filterwarnings('ignore')
import matplotlib
matplotlib.use('agg')
import numpy as np
import XFELDetAna.xfelpyanatools as xana
import XFELDetAna.xfelpycaltools as xcal
from cal_tools.enums import BadPixels
from cal_tools.tools import get_dir_creation_date
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from XFELDetAna.util import env
env.iprofile = cluster_profile
from XFELDetAna.detectors.jungfrau import readerPSI as jfreaderPSI
from XFELDetAna.detectors.jungfrau import reader as jfreader
from XFELDetAna.detectors.jungfrau.jf_chunk_reader import JFChunkReader
from XFELDetAna.detectors.jungfrau.util import non_empty_trains
# so constants relevant for the analysis
run_nums = [run_high, run_med, run_low] # run number for G0/HG0, G1, G2
cpuCores = 1
is_parallel = False
sensorSize = [1024, 512]
blockSize = [1024, 512]
xRange = [0, 0+sensorSize[0]]
yRange = [0, 0+sensorSize[1]]
gains = [0, 1, 2]
h5path = h5path.format(karabo_id, receiver_id)
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run_high)
print("Using {} as creation time".format(creation_time))
offset_abs_threshold = [offset_abs_threshold_low, offset_abs_threshold_high]
```
%% Cell type:code id: tags:
``` python
noiseCal = xcal.NoiseCalculator(sensorSize, nCells=memoryCells, cores=cpuCores, parallel=is_parallel, gains=gains,
blockSize=blockSize)
```
%% Cell type:code id: tags:
``` python
import h5py
for i_run, run in enumerate(run_nums):
if run is not None:
ped_dir = "{}/r{:04d}".format(in_folder, run)
fp_name = path_template.format(run, path_inset_control)
fp_path = '{}/{}'.format(ped_dir, fp_name)
if not manual_slow_data:
with h5py.File(fp_path.format(0), 'r') as f:
integration_time = int(f['/RUN/{}/DET/{}/exposureTime/value'.format(karabo_id, control_id)][()]*1e6)
bias_voltage = int(np.squeeze(f['/RUN/{}/DET/{}/vHighVoltage/value'.format(karabo_id, control_id)])[0])
print("Integration time is {} us".format(integration_time))
print("Bias voltage is {} V".format(bias_voltage))
fp_name = path_template.format(run, path_inset)
fp_path = '{}/{}'.format(ped_dir, fp_name)
filep_size = 500
myRange_P = range(0, sequences)
path = h5path
print("Reading data from {}".format(fp_path))
print("Run is: {}".format(run))
print("HDF5 path: {}".format(h5path))
reader = JFChunkReader(filename = fp_path, readFun = jfreader.readData, size = filep_size, chunkSize = chunkSize,
path = path, image_range=imageRange, pixels_x = sensorSize[0], pixels_y = sensorSize[1],
x_range = xRange, y_range = yRange, imagesPerChunk=chunkSize, filesRange = myRange_P,
memoryCells=memoryCells, blockSize=blockSize)
for data in reader.readChunks():
images = np.squeeze(data[0])
gainmaps = np.squeeze(data[1])
trainId = np.array(data[2])
idxs = non_empty_trains(trainId)
noiseCal.fill(data=images[..., idxs], gainMap=gainmaps[..., idxs])
else:
print("dark run for G{} is missing".format(i_run))
offset_map = noiseCal.getOffset()
noise_map = noiseCal.get()
```
%% Cell type:markdown id: tags:
## Offset and Noise Maps ##
Below offset and noise maps for the high ($g_0$) gain stage are shown, alongside the distribution of these values. One expects block-like structures mapping to the ASICs of the detector
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
%matplotlib inline
from XFELDetAna.plotting.histogram import histPlot
from XFELDetAna.plotting.heatmap import heatmapPlot
g_name = ['G0', 'G1', 'G2']
g_range = [(0, 8000), (8000, 16000), (8000, 16000)]
for g_idx in gains:
off = offset_map[:, :, 0, g_idx]
fo_0 = heatmapPlot(np.swapaxes(off, 0, 1),
y_label="Row",
x_label="Column",
lut_label="Pedestal {} [ADCu]".format(g_name[g_idx]),
vmin=g_range[g_idx][0],
vmax=g_range[g_idx][1])
fo_01, axo_01 = plt.subplots()
ho0, eo0 , xo0, stato0 = histPlot(axo_01, off,
bins=800,
range=g_range[g_idx],
facecolor='b',
histotype='stepfilled')
axo_01.tick_params(axis='both',which='major',labelsize=15)
axo_01.set_title('Module pedestal distribution', fontsize=15)
axo_01.set_xlabel('Pedestal {} [ADCu]'.format(g_name[g_idx]),fontsize=15)
axo_01.set_yscale('log')
noise = noise_map[:, :, 0, g_idx]
fn_0 = heatmapPlot(np.swapaxes(noise, 0, 1),
y_label="Row",
x_label="Column",
lut_label="RMS noise {} [ADCu]".format(g_name[g_idx]),
vmin=0,
vmax=50)
fn_01, axn_01 = plt.subplots()
hn0, en0 , xn0, statn0 = histPlot(axn_01, noise,
bins=100,
range=(0, 50),
facecolor='b',
histotype='stepfilled')
axn_01.tick_params(axis='both',which='major',labelsize=15)
axn_01.set_title('Module noise distribution', fontsize=15)
axn_01.set_xlabel('RMS noise {} [ADCu]'.format(g_name[g_idx]),fontsize=15)
axn_01.set_yscale('log')
plt.show()
```
%% Cell type:markdown id: tags:
## Bad Pixel Map ###
The bad pixel map is deduced by comparing offset and noise of each pixel ($v_i$) and each gain ($g$) against the median value for that gain stage:
$$
v_i > \mathrm{median}(v_{k,g}) + n \sigma_{v_{k,g}}
$$
or
$$
v_i < \mathrm{median}(v_{k,g}) - n \sigma_{v_{k,g}}
$$
Values are encoded in a 32 bit mask, where for the dark image deduced bad pixels the following non-zero entries are relevant:
%% Cell type:code id: tags:
``` python
def print_bp_entry(bp):
print("{:<30s} {:032b}".format(bp.name, bp.value))
print_bp_entry(BadPixels.OFFSET_OUT_OF_THRESHOLD)
print_bp_entry(BadPixels.NOISE_OUT_OF_THRESHOLD)
print_bp_entry(BadPixels.OFFSET_NOISE_EVAL_ERROR)
```
%% Cell type:code id: tags:
``` python
bad_pixels_map = np.zeros(noise_map.shape, np.uint32)
def eval_bpidx(d):
mdn = np.nanmedian(d, axis=(0, 1, 2))[None, None, None, :]
std = np.nanstd(d, axis=(0, 1, 2))[None, None, None, :]
idx = (d > badpixel_threshold_sigma*std+mdn) | (d < (-badpixel_threshold_sigma)*std+mdn)
return idx
offset_abs_threshold = np.array(offset_abs_threshold)
bad_pixels_map[eval_bpidx(offset_map)] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value
bad_pixels_map[~np.isfinite(offset_map)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
bad_pixels_map[eval_bpidx(noise_map)] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
bad_pixels_map[~np.isfinite(noise_map)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
bad_pixels_map[(offset_map < offset_abs_threshold[0][None, None, None, :]) | (offset_map > offset_abs_threshold[1][None, None, None, :])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value
for g_idx in gains:
bad_pixels = bad_pixels_map[:, :, 0, g_idx]
fn_0 = heatmapPlot(np.swapaxes(bad_pixels, 0, 1),
y_label="Row",
x_label="Column",
lut_label="Badpixels {} [ADCu]".format(g_name[g_idx]),
vmin=0)
```
%% Cell type:code id: tags:
``` python
## offset
metadata = ConstantMetaData()
offset = Constants.jungfrau.Offset()
offset.data = np.moveaxis(offset_map, 0, 1)
metadata.calibration_constant = offset
# set the operating condition
condition = Conditions.Dark.jungfrau(memory_cells=1, bias_voltage=bias_voltage,
integration_time=integration_time)
device = getattr(Detectors, db_module)
metadata.detector_condition = condition
# specify the a version for this constant
if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device)
metadata.retrieve(cal_db_interface)
else:
metadata.calibration_constant_version = Versions.Timespan(device=device,
start=creation_time)
metadata.send(cal_db_interface)
## noise
metadata = ConstantMetaData()
noise = Constants.jungfrau.Noise()
noise.data = np.moveaxis(noise_map, 0, 1)
metadata.calibration_constant = noise
# set the operating condition
condition = Conditions.Dark.jungfrau(memory_cells=1, bias_voltage=bias_voltage,
integration_time=integration_time)
metadata.detector_condition = condition
# specify the a version for this constant
if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device)
metadata.retrieve(cal_db_interface)
else:
metadata.calibration_constant_version = Versions.Timespan(device=device,
start=creation_time)
metadata.send(cal_db_interface)
## bad pixels
metadata = ConstantMetaData()
bpix = Constants.jungfrau.BadPixelsDark()
bpix.data = np.moveaxis(bad_pixels_map, 0, 1)
metadata.calibration_constant = bpix
# set the operating condition
condition = Conditions.Dark.jungfrau(memory_cells=1, bias_voltage=bias_voltage,
integration_time=integration_time)
metadata.detector_condition = condition
# specify the a version for this constant
if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device)
metadata.retrieve(cal_db_interface)
else:
metadata.calibration_constant_version = Versions.Timespan(device=device,
start=creation_time)
metadata.send(cal_db_interface)
```
%% Cell type:markdown id: tags:
## Verification ##
For verification we correct the input data with the offset deduced from it. We expect a distribution centered around zero.
%% Cell type:code id: tags:
``` python
fd_path = fp_path
d_path = path
filed_size = 500
myRange_D = range(0, 2)
imageRange = [0, 500]
chunkSize = 25
reader = JFChunkReader(filename = fd_path, readFun = jfreader.readData, size = filed_size, chunkSize = chunkSize,
path = d_path, image_range=imageRange, pixels_x = sensorSize[0], pixels_y = sensorSize[1],
x_range = xRange, y_range = yRange, imagesPerChunk=chunkSize, filesRange = myRange_D,
memoryCells=memoryCells, blockSize=blockSize)
offsetCorr = xcal.OffsetCorrection(shape=sensorSize, offsetMap=offset_map, nCells=memoryCells,
cores=cpuCores, parallel=is_parallel, gains=gains, blockSize=blockSize)
pixel_data = []
for data_raw in reader.readChunks():
images = np.array(data_raw[0]).astype(np.float32)
gainmaps = np.array(data_raw[1])
trainID = np.array(data_raw[2])
idxs = non_empty_trains(trainId)
data_out = offsetCorr.correct(data=images[..., idxs], gainMap=gainmaps[..., idxs])
pixel_data.append(data_out)
pixel_data = np.array(pixel_data)
```
%% Cell type:code id: tags:
``` python
is_log = True
h_range = (-500, 500)
h_bins = 300
f_sp, ax_sp = plt.subplots()
h, e , x, stat = histPlot(ax_sp, pixel_data.flatten(),
bins=h_bins,
range=h_range,
facecolor='b',
log=is_log,
histotype='stepfilled')
ax_sp.tick_params(axis='both',which='major',labelsize=15)
ax_sp.set_xlabel('Energy [ADCu]', fontsize=15)
plt.show()
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Legal Notice | Contacts