Skip to content
Snippets Groups Projects
Commit 820796d2 authored by Mikhail Karnevskiy's avatar Mikhail Karnevskiy
Browse files

Fix title

parent 72c49314
No related branches found
No related tags found
1 merge request!102Feat: integration of epix10K
Expand all Hide whitespace changes
Inline Side-by-side
notebooks/ePix/Characterize_Darks_ePix10K_NBC.ipynb
+ 14
− 92
View file @ 820796d2
This diff is collapsed.
notebooks/ePix/Correction_ePix10K_NBC.ipynb
+ 1
− 1
View file @ 820796d2
%% Cell type:markdown id: tags:
# ePIX10K Data Correction ##
# ePIX10K Data Correction #
Authors: M. Karnevskiy, S. Hauf, Version 1.0
The following notebook provides Offset correction of images acquired with the ePix10K detector.
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/HED/201922/p002550/raw" # input folder, required
out_folder = '/gpfs/exfel/data/scratch/karnem/test/' # output folder, required
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data
run = 55 # which run to read data from, required
h5path = '/INSTRUMENT/{}/DET/RECEIVER:daqOutput/data/image' # path in the HDF5 file to images
h5path_t = '/INSTRUMENT/{}/DET/RECEIVER:daqOutput/data/backTemp' # path to find temperature at
h5path_cntrl = '/CONTROL/{}/DET' # path to control data
path_inset = "EPIX03" # file inset for image data
instance = "HED_IA1_EPIX10K-1" # karabo instance
cluster_profile = "noDB" # ipcluster profile to use
cpuCores = 4 # Specifies the number of running cpu cores
cal_db_interface = "tcp://max-exfl016:8015#8025" # calibration DB interface to use
cal_db_timeout = 30000000 # timeout on caldb requests
sequences = [-1] # sequences to correct, set to -1 for all, range allowed
chunk_size_idim = 1 # H5 chunking size of output data
overwrite = True # overwrite output folder
limit_images = 0 # process only first N images, 0 - process all
use_dir_creation_date = True # date constants injected before directory creation time
db_module = "ePix10K_M43" # module id in the database
sequences_per_node = 1 # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel
bias_voltage = 200 # bias voltage
in_vacuum = False # detector operated in vacuum
fix_temperature = 290. # fix temperature to this value
split_evt_primary_threshold = 7. # primary threshold for split event correction
split_evt_secondary_threshold = 5. # secondary threshold for split event correction
split_evt_mip_threshold = 1000. # minimum ionizing particle threshold
def balance_sequences(in_folder, run, sequences, sequences_per_node):
import glob
import re
import numpy as np
if sequences[0] == -1:
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)
else:
seq_nums = set(sequences)
nsplits = len(seq_nums)//sequences_per_node+1
while nsplits > 8:
sequences_per_node += 1
nsplits = len(seq_nums)//sequences_per_node+1
print("Changed to {} sequences per node to have a maximum of 8 concurrent jobs".format(sequences_per_node))
return [l.tolist() for l in np.array_split(list(seq_nums), nsplits) if l.size > 0]
```
%% Cell type:code id: tags:
``` python
import XFELDetAna.xfelprofiler as xprof
profiler = xprof.Profiler()
profiler.disable()
from XFELDetAna.util import env
env.iprofile = cluster_profile
import warnings
warnings.filterwarnings('ignore')
from XFELDetAna import xfelpycaltools as xcal
from XFELDetAna import xfelpyanatools as xana
from XFELDetAna.plotting.util import prettyPlotting
prettyPlotting=True
from XFELDetAna.xfelreaders import ChunkReader
from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5
import numpy as np
import h5py
import time
import copy
import os
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from iCalibrationDB.detectors import DetectorTypes
from cal_tools.tools import get_dir_creation_date, get_constant_from_db
%matplotlib inline
if sequences[0] == -1:
sequences = None
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("Interface to Cal DB: ", cal_db_interface)
h5path = h5path.format(instance)
h5path_t = h5path_t.format(instance)
h5path_cntrl = h5path_cntrl.format(instance)
```
%% Cell type:code id: tags:
``` python
x = 356 # rows of the ePix10K
y = 384 # columns of the ePix10K
ped_dir = "{}/r{:04d}".format(in_folder, run)
out_folder = "{}/r{:04d}".format(out_folder, run)
fp_name = path_template.format(run, path_inset)
fp_path = '{}/{}'.format(ped_dir, fp_name)
print("Reading from: {}".format(fp_path))
print("Run is: {}".format(run))
print("HDF5 path: {}".format(h5path))
print("Data is output to: {}".format(out_folder))
import datetime
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run)
if creation_time:
print("Using {} as creation time".format(creation_time.isoformat()))
```
%% Cell type:code id: tags:
``` python
sensorSize = [x, y]
chunkSize = 100 #Number of images to read per chunk
blockSize = [sensorSize[0]//2, sensorSize[1]//2] # Sensor area will be analysed according to blocksize
xcal.defaultBlockSize = blockSize
memoryCells = 1 # ePIX10K has no memory cell
run_parallel = True
filename = fp_path.format(sequences[0] if sequences else 0)
with h5py.File(filename, 'r') as f:
integration_time = int(f['{}/CONTROL/expTime/value'.format(h5path_cntrl)][0])
gain_setting = int(f['{}/CONTROL/encodedGainSetting/value'.format(h5path_cntrl)][0])
temperature = np.mean(f[h5path_t])/100.
temperature_k = temperature + 273.15
if fix_temperature != 0:
temperature_k = fix_temperature
print("Temperature is fixed!")
print("Bias voltage is {} V".format(bias_voltage))
print("Gain setting {}".format(gain_setting))
print("Detector integration time is set to {}".format(integration_time))
print("Mean temperature was {:0.2f} °C / {:0.2f} K at beginning of run".format(temperature, temperature_k))
print("Operated in vacuum: {} ".format(in_vacuum))
if not os.path.exists(out_folder):
os.makedirs(out_folder)
elif not overwrite:
raise AttributeError("Output path exists! Exiting")
```
%% 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
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))
if len(table):
md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=["#", "file"])))
```
%% Cell type:markdown id: tags:
As a first step, dark maps have to be loaded.
%% Cell type:code id: tags:
``` python
dclass = "ePix10K"
const_name = "Offset"
offsetMap = None
temp_limits = 5.
# offset
det = getattr(Detectors, db_module)
dconstants = getattr(Constants, dclass)
dcond = Conditions.Dark
condition = getattr(dcond, dclass)(bias_voltage=bias_voltage,
integration_time=integration_time,
temperature=temperature_k,
in_vacuum=in_vacuum,
gain_setting=gain_setting)
for parm in condition.parameters:
if parm.name == "Sensor Temperature":
parm.lower_deviation = temp_limits
parm.upper_deviation = temp_limits
offsetMap = get_constant_from_db(det,
getattr(dconstants, const_name)(),
condition,
None,
cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout)
condition = getattr(dcond, dclass)(bias_voltage=bias_voltage,
integration_time=integration_time,
temperature=temperature_k,
in_vacuum=in_vacuum,
gain_setting=gain_setting)
const_name = "Noise"
noiseMap = get_constant_from_db(det,
getattr(dconstants, const_name)(),
condition,
None,
cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout)
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
#************************Calculators************************#
offsetCorrection = xcal.OffsetCorrection(sensorSize,
offsetMap,
nCells = memoryCells,
cores=cpuCores, gains=None,
blockSize=blockSize,
parallel=run_parallel)
```
%% Cell type:code id: tags:
``` python
#*****************Histogram Calculators******************#
histCalOffsetCor = xcal.HistogramCalculator(sensorSize,
bins=1050,
range=[-50, 1000], parallel=run_parallel,
nCells=memoryCells,
cores=cpuCores,
blockSize=blockSize)
```
%% Cell type:markdown id: tags:
Applying corrections
%% Cell type:code id: tags:
``` python
histCalOffsetCor.debug()
offsetCorrection.debug()
```
%% Cell type:code id: tags:
``` python
#************************Calculators************************#
commonModeBlockSize = [x//2, y//2]
commonModeAxisR = 'row'
cmCorrection = xcal.CommonModeCorrection([x, y],
commonModeBlockSize,
commonModeAxisR,
nCells = memoryCells,
noiseMap = noiseMap,
runParallel=True,
stats=True)
patternClassifier = xcal.PatternClassifier([x, y],
noiseMap,
split_evt_primary_threshold,
split_evt_secondary_threshold,
split_evt_mip_threshold,
tagFirstSingles = 0,
nCells=memoryCells,
cores=cpuCores,
allowElongated = False,
blockSize=[x, y],
runParallel=True)
histCalCMCor = xcal.HistogramCalculator(sensorSize,
bins=1050,
range=[-50, 1000], parallel=run_parallel,
nCells=memoryCells,
cores=cpuCores,
blockSize=blockSize)
histCalSECor = xcal.HistogramCalculator(sensorSize,
bins=1050,
range=[-50, 1000], parallel=run_parallel,
nCells=memoryCells,
cores=cpuCores,
blockSize=blockSize)
```
%% Cell type:code id: tags:
``` python
cmCorrection.debug()
patternClassifier.debug()
histCalCMCor.debug()
histCalSECor.debug()
```
%% 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 `correct EPIX10K`
"""
if h5base.startswith("/"):
h5base = h5base[1:]
dont_copy = ['pixels']
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
for k, f in enumerate(file_list):
with h5py.File(f, 'r') as infile:
out_fileb = "{}/{}".format(out_folder, f.split("/")[-1])
out_file = out_fileb.replace("RAW", "CORR")
#out_filed = out_fileb.replace("RAW", "CORR-SC")
data = None
with h5py.File(out_file, "w") as ofile:
try:
copy_and_sanitize_non_cal_data(infile, ofile, h5path)
data = infile[h5path+"/pixels"][()]
data = np.compress(np.any(data>0, axis=(1,2)), data, axis=0)
if limit_images > 0:
data = data[:limit_images,...]
oshape = data.shape
data = np.moveaxis(data, 0, 2)
ddset = ofile.create_dataset(h5path+"/pixels",
oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.float32)
data = offsetCorrection.correct(data.astype(np.float32)) #correct for the offset
histCalOffsetCor.fill(data)
ddset[...] = np.moveaxis(data, 2, 0)
except Exception as e:
print("Couldn't calibrate data in {}: {}".format(f, e))
if False:
with h5py.File(out_file, "a") as ofiled:
try:
#copy_and_sanitize_non_cal_data(infile, ofiled, h5path)
ddsetcm = ofiled.create_dataset(h5path+"/pixels_cm",
oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.float32)
ddsetc = ofiled.create_dataset(h5path+"/pixels_classified",
oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.float32, compression="gzip")
ddsetp = ofiled.create_dataset(h5path+"/patterns",
oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.int32, compression="gzip")
data = cmCorrection.correct(data) #correct for the row common mode
histCalCMCor.fill(data)
ddsetcm[...] = np.moveaxis(data, 2, 0)
data, patterns = patternClassifier.classify(data)
data[data < (split_evt_primary_threshold*noiseMap)] = 0
ddsetc[...] = np.moveaxis(data, 2, 0)
ddsetp[...] = np.moveaxis(patterns, 2, 0)
data[patterns != 100] = np.nan
histCalSECor.fill(data)
except Exception as e:
print("Couldn't calibrate data in {}: {}".format(f, e))
```
%% Cell type:code id: tags:
``` python
ho,eo,co,so = histCalOffsetCor.get()
d = [{'x': co,
'y': ho,
'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'Offset corr.'
},
]
if False:
ho,eo,co,so = histCalCMCor.get()
d.append({'x': co,
'y': ho,
'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'CM corr.'
})
ho,eo,co,so = histCalSECor.get()
d.append({'x': co,
'y': ho,
'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'Single split events'
})
fig = xana.simplePlot(d, aspect=1, x_label='Energy(ADU)',
y_label='Number of occurrences', figsize='2col',
y_log=True, x_range=(-50,500),
legend='top-center-frame-2col')
```
%% Cell type:markdown id: tags:
## Mean Image of last Sequence ##
%% Cell type:code id: tags:
``` python
fig = xana.heatmapPlot(np.nanmedian(data, axis=2),
x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)',
x_range=(0,y),
y_range=(0,x), vmin=-50, vmax=50)
```
%% Cell type:markdown id: tags:
## Single Shot of last Sequence ##
%% Cell type:code id: tags:
``` python
fig = xana.heatmapPlot(data[...,0],
x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)',
x_range=(0,y),
y_range=(0,x), vmin=-50, vmax=50)
```
%% Cell type:code id: tags:
``` python
```
......
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