Skip to content
Snippets Groups Projects
Commit bbc35c19 authored by Michele Cascella's avatar Michele Cascella Committed by Karim Ahmed
Browse files

add common mode correction to epix correction notebook + WIP: patter classification.

parent 5d46f8e6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
notebooks/ePix100/Correction_ePix100_NBC.ipynb
+ 116
− 70
View file @ bbc35c19
%% Cell type:markdown id: tags:
# ePIX Data Correction ##
# ePIX Data Correction
Authors: Q. Tian S. Hauf, Version 1.0
Authors: Q. Tian S. Hauf M. Cascella, Version 1.0
The following notebook provides Offset correction of images acquired with the ePix100 detector.
%% Cell type:code id: tags:
``` python
cluster_profile = "noDB" # ipcluster profile to use
in_folder = "/gpfs/exfel/exp/CALLAB/202031/p900113/raw" # input folder, required
out_folder = "" # output folder, required
sequences = [-1] # sequences to correct, set to -1 for all, range allowed
run = 9988 # which run to read data from, required
karabo_id = "MID_EXP_EPIX-1" # karabo karabo_id
karabo_da = "EPIX01" # data aggregators
db_module = "ePix100_M15" # module id in the database
receiver_id = "RECEIVER" # inset for receiver devices
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data
h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data/image' # path in the HDF5 file to images
h5path_t = '/INSTRUMENT/{}/DET/{}:daqOutput/data/backTemp' # path to find temperature at
h5path_cntrl = '/CONTROL/{}/DET' # path to control data
use_dir_creation_date = True # date constants injected before directory creation time
cal_db_interface = "tcp://max-exfl016:8015#8025" # calibration DB interface to use
cal_db_timeout = 300000 # timeout on caldb requests
cpuCores = 4 # Specifies the number of running cpu cores
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
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
gain_photon_energy = 9.0 # Photon energy used for gain calibration
photon_energy = 8.0 # Photon energy to calibrate in number of photons, 0 for calibration in keV
relative_gain = False # Apply relative gain correction.
common_mode = True # Apply common mode correction.
cm_min_frac = 0.25 # No CM correction is performed if after masking the ratio of good pixels falls below this
cm_noise_sigma = 5. # CM correction noise standard deviation
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, karabo_da):
from xfel_calibrate.calibrate import balance_sequences as bs
return bs(in_folder, run, sequences, sequences_per_node, karabo_da)
```
%% Cell type:code id: tags:
``` python
import tabulate
import warnings
import h5py
import numpy as np
from IPython.display import Latex, display
from pathlib import Path
import XFELDetAna.xfelprofiler as xprof
from XFELDetAna import xfelpyanatools as xana
from XFELDetAna import xfelpycaltools as xcal
from XFELDetAna.plotting.util import prettyPlotting
from XFELDetAna.util import env
from cal_tools.tools import (
get_constant_from_db,
get_dir_creation_date,
)
from iCalibrationDB import (
Conditions,
Constants,
)
warnings.filterwarnings('ignore')
prettyPlotting = True
profiler = xprof.Profiler()
profiler.disable()
env.iprofile = cluster_profile
%matplotlib inline
```
%% Cell type:code id: tags:
``` python
# TODO: expose to first cell after fixing common mode correction.
common_mode = False # Apply common mode correction.
# TODO: expose to first cell after fixing clustering.
pattern_classification = False # do clustering.
h5path = h5path.format(karabo_id, receiver_id)
h5path_t = h5path_t.format(karabo_id, receiver_id)
h5path_cntrl = h5path_cntrl.format(karabo_id)
plot_unit = 'ADU'
if relative_gain:
plot_unit = 'keV'
if photon_energy > 0:
plot_unit = '$\gamma$'
```
%% Cell type:code id: tags:
``` python
x = 708 # rows of the ePix100
y = 768 # columns of the ePix100
in_folder = Path(in_folder)
ped_dir = in_folder / f"r{run:04d}"
fp_name = path_template.format(run, karabo_da)
print(f"Reading from: {ped_dir / fp_name}")
print(f"Run is: {run}")
print(f"HDF5 path: {h5path}")
print(f"Data is output to: {out_folder}")
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run)
if creation_time:
print(f"Using {creation_time.isoformat()} as creation time")
```
%% Cell type:code id: tags:
``` python
sensorSize = [x, y]
chunkSize = 100 # Number of images to read per chunk
# Sensor area will be analysed according to blocksize
blockSize = [sensorSize[0]//2, sensorSize[1]//2]
xcal.defaultBlockSize = blockSize
memoryCells = 1 # ePIX has no memory cells
run_parallel = True
# Read slow data from the first available sequence file.
filename = ped_dir / fp_name.format(sequences[0] if sequences[0] != -1 else 0)
filename = ped_dir / fp_name.format(
sequences[0] if sequences[0] != -1 else 0)
with h5py.File(filename, 'r') as f:
integration_time = int(f[f"{h5path_cntrl}/CONTROL/expTime/value"][0])
integration_time = int(
f[f"{h5path_cntrl}/CONTROL/expTime/value"][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(f"Bias voltage is {bias_voltage} V")
print(f"Detector integration time is set to {integration_time}")
print(
f"Mean temperature was {temperature:0.2f} °C "
f"/ {temperature_k:0.2f} K at beginning of run"
)
print(f"Operated in vacuum: {in_vacuum} ")
out_folder = Path(out_folder)
if out_folder.is_dir() and not overwrite:
raise AttributeError("Output path exists! Exiting")
out_folder.mkdir(parents=True, exist_ok=True)
```
%% Cell type:code id: tags:
``` python
# Glob the right *.h5 fast data files.
seq_files = sorted(ped_dir.glob(f"*{karabo_da}*.h5"))
# If a set of sequences requested to correct,
# adapt seq_files list.
if sequences != [-1]:
seq_files = [f for f in seq_files
if any(f.match(f"*-S{s:05d}.h5") for s in sequences)]
if any(
f.match(f"*-S{s:05d}.h5") for s in sequences)]
print(f"Processing a total of {len(seq_files)} sequence files")
```
%% Cell type:code id: tags:
``` python
# Table of sequence files to process
table = [(k, f) for k, f in enumerate(seq_files)]
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
temp_limits = 5.
cond_dict = {
"bias_voltage": bias_voltage,
"integration_time": integration_time,
"temperature": temperature_k,
"in_vacuum": in_vacuum,
}
dark_condition = Conditions.Dark.ePix100(**cond_dict)
# update conditions with illuminated conditins.
cond_dict.update({
"photon_energy": gain_photon_energy
})
illum_condition = Conditions.Illuminated.ePix100(**cond_dict)
const_cond = {
"Offset": dark_condition,
"Noise": dark_condition,
"RelativeGain": illum_condition,
}
const_data = dict()
for cname, condition in const_cond.items():
if cname == "RelativeGain" and not relative_gain:
continue
# TODO: Fix this logic.
for parm in condition.parameters:
if parm.name == "Sensor Temperature":
parm.lower_deviation = temp_limits
parm.upper_deviation = temp_limits
const_data[cname] = get_constant_from_db(
karabo_id=karabo_id,
karabo_da=karabo_da,
constant=getattr(Constants.ePix100, cname)(),
condition=condition,
empty_constant=None,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
print_once=2,
timeout=cal_db_timeout
)
if relative_gain and const_data["RelativeGain"] is None:
print(
"WARNING: RelativeGain map is requested, but not found./n"
"No gain correction will be applied"
)
relative_gain = False
plot_unit = 'ADU'
```
%% Cell type:code id: tags:
``` python
# ************************Calculators************************ #
# ************************Calculators******************** #
offsetCorrection = xcal.OffsetCorrection(
sensorSize,
const_data["Offset"],
nCells=memoryCells,
cores=cpuCores,
gains=None,
blockSize=blockSize,
parallel=run_parallel
)
if relative_gain:
gainCorrection = xcal.RelativeGainCorrection(
sensorSize,
1./const_data["RelativeGain"][..., None],
nCells=memoryCells,
parallel=run_parallel,
cores=cpuCores,
blockSize=blockSize,
gains=None,
)
```
%% 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()
if relative_gain:
gainCorrection.debug()
```
%% Cell type:code id: tags:
``` python
# ************************Calculators************************ #
# ************************Calculators******************** #
if common_mode:
commonModeBlockSize = [x//2, y//2]
commonModeAxisR = 'row'
cmCorrection = xcal.CommonModeCorrection(
[x, y],
commonModeBlockSize,
commonModeAxisR,
stats = True
cmCorrectionB = xcal.CommonModeCorrection(
shape=sensorSize,
blockSize=commonModeBlockSize,
orientation='block',
nCells=memoryCells,
noiseMap=const_data["Noise"],
noiseMap=const_data['Noise'],
runParallel=run_parallel,
stats=True,
stats=stats,
minFrac=cm_min_frac,
noiseSigma=cm_noise_sigma,
)
cmCorrectionR = xcal.CommonModeCorrection(
shape=sensorSize,
blockSize=commonModeBlockSize,
orientation='row',
nCells=memoryCells,
noiseMap=const_data['Noise'],
runParallel=run_parallel,
stats=stats,
minFrac=cm_min_frac,
noiseSigma=cm_noise_sigma,
)
cmCorrectionC = xcal.CommonModeCorrection(
shape=sensorSize,
blockSize=commonModeBlockSize,
orientation='col',
nCells=memoryCells,
noiseMap=const_data['Noise'],
runParallel=run_parallel,
stats=stats,
minFrac=cm_min_frac,
noiseSigma=cm_noise_sigma,
)
histCalCMCor = xcal.HistogramCalculator(
sensorSize,
bins=1050,
range=[-50, 1000],
parallel=run_parallel,
nCells=memoryCells,
cores=cpuCores,
blockSize=blockSize,
)
```
patternClassifier = xcal.PatternClassifier(
[x, y],
const_data["Noise"],
split_evt_primary_threshold,
split_evt_secondary_threshold,
split_evt_mip_threshold,
tagFirstSingles=0,
nCells=memoryCells,
cores=cpuCores,
allowElongated=False,
blockSize=[x, y],
runParallel=run_parallel,
)
%% Cell type:code id: tags:
histCalSECor = xcal.HistogramCalculator(
sensorSize,
bins=1050,
range=[-50, 1000],
parallel=run_parallel,
nCells=memoryCells,
cores=cpuCores,
blockSize=blockSize,
)
``` python
if pattern_classification:
patternClassifier = xcal.PatternClassifier(
[x, y],
const_data["Noise"],
split_evt_primary_threshold,
split_evt_secondary_threshold,
split_evt_mip_threshold,
tagFirstSingles=0,
nCells=memoryCells,
cores=cpuCores,
allowElongated=False,
blockSize=[x, y],
runParallel=run_parallel,
)
histCalSECor = xcal.HistogramCalculator(
sensorSize,
bins=1050,
range=[-50, 1000],
parallel=run_parallel,
nCells=memoryCells,
cores=cpuCores,
blockSize=blockSize,
)
```
%% Cell type:code id: tags:
``` python
if common_mode:
cmCorrection.debug()
cmCorrectionB.debug()
cmCorrectionR.debug()
cmCorrectionC.debug()
histCalCMCor.debug()
patternClassifier.debug()
histCalSECor.debug()
if pattern_classification:
patternClassifier.debug()
histCalSECor.debug()
```
%% Cell type:code id: tags:
``` python
def copy_and_sanitize_non_cal_data(
infile: h5py,
outfile: h5py,
h5base: str
):
""" Copy and sanitize data in `infile`
that is not touched by `correctEPIX`. """
if h5base.startswith("/"):
h5base = h5base[1:]
dont_copy = [h5base+"/pixels"]
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 f in seq_files:
data = None
out_file = out_folder / f.name.replace("RAW", "CORR")
with h5py.File(f, "r") as infile, h5py.File(out_file, "w") as ofile:
with h5py.File(f, "r") as infile, h5py.File(out_file, "w") as ofile: # noqa
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)
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)
# Offset correction.
data = offsetCorrection.correct(data.astype(np.float32))
# Common Mode correction.
if common_mode:
# Block CM
data = cmCorrectionB.correct(data)
# Row CM
data = cmCorrectionR.correct(data)
# COL CM
data = cmCorrectionC.correct(data)
histCalCMCor.fill(data)
# relative gain correction.
if relative_gain:
data = gainCorrection.correct(data.astype(np.float32))
if photon_energy > 0:
data /= photon_energy
histCalOffsetCor.fill(data)
ddset[...] = np.moveaxis(data, 2, 0)
# Common mode correction
# TODO: Fix conflict between common mode and relative gain.
"""The gain correction is currently applying an absolute correction
(not a relative correction as the implied by the name);
it changes the scale (the unit of measurement) of the data from ADU
to either keV or n_of_photons. But the common mode correction
relies on comparing data with the noise map, which is still in ADU.
The best solution is probably to do a relative gain correction first
(correct) and apply the global absolute gain to the data at the end,
after common mode and clustering.
"""The gain correction is currently applying
an absolute correction (not a relative correction
as the implied by the name);
it changes the scale (the unit of measurement)
of the data from ADU to either keV or n_of_photons.
But the pattern classification relies on comparing
data with the noise map, which is still in ADU.
The best solution is to do a relative gain
correction first and apply the global absolute
gain to the data at the end, after clustering.
"""
if common_mode:
ddsetcm = ofile.create_dataset(
h5path+"/pixels_cm",
oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.float32)
# TODO: Fix conflict between pattern classification
# and gain corr.
if pattern_classification:
ddsetc = ofile.create_dataset(
h5path+"/pixels_classified",
oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.float32, compression="gzip")
ddsetp = ofile.create_dataset(
h5path+"/patterns",
oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.int32, compression="gzip")
# row common mode correction.
data = cmCorrection.correct(data)
histCalCMCor.fill(data)
ddsetcm[...] = np.moveaxis(data, 2, 0)
data, patterns = patternClassifier.classify(data)
data[data < (split_evt_primary_threshold*const_data["Noise"])] = 0 # noqa
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(f"ERROR applying common mode correction for {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 common_mode:
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.'
})
if pattern_classification:
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=f'Energy({plot_unit})',
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=f'Signal ({plot_unit})',
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=f'Signal ({plot_unit})',
x_range=(0, y),
y_range=(0, x),
vmin=-50, vmax=50)
```
......
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