Skip to content
Snippets Groups Projects
Commit cb46bba3 authored by Karim Ahmed's avatar Karim Ahmed
Browse files

remove unrelated changes

parent b0da6ce9
No related branches found
No related tags found
1 merge request!783[JUNGFRAU][CORRECT] Fix/store roi only if defined
Hide whitespace changes
Inline Side-by-side
notebooks/ePix100/Correction_ePix100_NBC.ipynb
+ 63
− 97
View file @ cb46bba3
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# ePix100 Data Correction # ePix100 Data Correction
Author: European XFEL Detector Group, Version: 2.0 Author: European XFEL Detector Group, Version: 2.0
The following notebook provides data correction of images acquired with the ePix100 detector. The following notebook provides data correction of images acquired with the ePix100 detector.
The sequence of correction applied are: The sequence of correction applied are:
Offset --> Common Mode Noise --> Relative Gain --> Charge Sharing --> Absolute Gain. Offset --> Common Mode Noise --> Relative Gain --> Charge Sharing --> Absolute Gain.
Offset, common mode and gain corrected data is saved to /data/image/pixels in the CORR files. Offset, common mode and gain corrected data is saved to /data/image/pixels in the CORR files.
If pattern classification is applied (charge sharing correction), this data will be saved to /data/image/pixels_classified, while the corresponding patterns will be saved to /data/image/patterns in the CORR files. If pattern classification is applied (charge sharing correction), this data will be saved to /data/image/pixels_classified, while the corresponding patterns will be saved to /data/image/patterns in the CORR files.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
in_folder = "/gpfs/exfel/exp/HED/202202/p003121/raw" # input folder, required in_folder = "/gpfs/exfel/exp/HED/202202/p003121/raw" # input folder, required
out_folder = "" # output folder, required out_folder = "" # output folder, required
metadata_folder = "" # Directory containing calibration_metadata.yml when run by xfel-calibrate metadata_folder = "" # Directory containing calibration_metadata.yml when run by xfel-calibrate
sequences = [-1] # sequences to correct, set to -1 for all, range allowed sequences = [-1] # sequences to correct, set to -1 for all, range allowed
sequences_per_node = 1 # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel sequences_per_node = 1 # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel
run = 156 # which run to read data from, required run = 156 # which run to read data from, required
# Parameters for accessing the raw data. # Parameters for accessing the raw data.
karabo_id = "HED_IA1_EPX100-1" # karabo karabo_id karabo_id = "HED_IA1_EPX100-1" # karabo karabo_id
karabo_da = "EPIX01" # data aggregators karabo_da = "EPIX01" # data aggregators
db_module = "" # module id in the database db_module = "" # module id in the database
receiver_template = "RECEIVER" # detector receiver template for accessing raw data files receiver_template = "RECEIVER" # detector receiver template for accessing raw data files
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data
instrument_source_template = '{}/DET/{}:daqOutput' # instrument detector data source in h5files instrument_source_template = '{}/DET/{}:daqOutput' # instrument detector data source in h5files
# Parameters affecting writing corrected data. # Parameters affecting writing corrected data.
chunk_size_idim = 1 # H5 chunking size of output data chunk_size_idim = 1 # H5 chunking size of output data
limit_trains = 0 # Process only first N images, 0 - process all.
# Only for testing
limit_images = 0 # ONLY FOR TESTING. process only first N images, 0 - process all.
# Parameters for the calibration database. # Parameters for the calibration database.
cal_db_interface = "tcp://max-exfl016:8015#8025" # calibration DB interface to use cal_db_interface = "tcp://max-exfl016:8015#8025" # calibration DB interface to use
cal_db_timeout = 300000 # timeout on caldb requests cal_db_timeout = 300000 # timeout on caldb requests
creation_time = "" # The timestamp to use with Calibration DBe. Required Format: "YYYY-MM-DD hh:mm:ss" e.g. 2019-07-04 11:02:41 creation_time = "" # The timestamp to use with Calibration DBe. Required Format: "YYYY-MM-DD hh:mm:ss" e.g. 2019-07-04 11:02:41
# Conditions for retrieving calibration constants. # Conditions for retrieving calibration constants.
bias_voltage = 200 # bias voltage bias_voltage = 200 # bias voltage
in_vacuum = False # detector operated in vacuum in_vacuum = False # detector operated in vacuum
integration_time = -1 # Detector integration time, Default value -1 to use the value from the slow data. integration_time = -1 # Detector integration time, Default value -1 to use the value from the slow data.
fix_temperature = -1 # fixed temperature value in Kelvin, Default value -1 to use the value from files. fix_temperature = -1 # fixed temperature value in Kelvin, Default value -1 to use the value from files.
gain_photon_energy = 8.048 # Photon energy used for gain calibration gain_photon_energy = 8.048 # Photon energy used for gain calibration
photon_energy = 0. # Photon energy to calibrate in number of photons, 0 for calibration in keV photon_energy = 0. # Photon energy to calibrate in number of photons, 0 for calibration in keV
# Flags to select type of applied corrections. # Flags to select type of applied corrections.
pattern_classification = True # do clustering. pattern_classification = True # do clustering.
relative_gain = True # Apply relative gain correction. relative_gain = True # Apply relative gain correction.
absolute_gain = True # Apply absolute gain correction (implies relative gain). absolute_gain = True # Apply absolute gain correction (implies relative gain).
common_mode = True # Apply common mode correction. common_mode = True # Apply common mode correction.
# Parameters affecting applied correction. # Parameters affecting applied correction.
cm_min_frac = 0.25 # No CM correction is performed if after masking the ratio of good pixels falls below this 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 cm_noise_sigma = 5. # CM correction noise standard deviation
split_evt_primary_threshold = 7. # primary threshold for split event correction split_evt_primary_threshold = 7. # primary threshold for split event correction
split_evt_secondary_threshold = 5. # secondary 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 split_evt_mip_threshold = 1000. # minimum ionizing particle threshold
def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da): def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da):
from xfel_calibrate.calibrate import balance_sequences as bs from xfel_calibrate.calibrate import balance_sequences as bs
return bs(in_folder, run, sequences, sequences_per_node, karabo_da) return bs(in_folder, run, sequences, sequences_per_node, karabo_da)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import tabulate import tabulate
import warnings import warnings
from logging import warning
from sys import exit
import h5py import h5py
import pasha as psh import pasha as psh
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from IPython.display import Latex, display from IPython.display import Latex, display
from extra_data import RunDirectory, H5File from extra_data import RunDirectory, H5File
from pathlib import Path from pathlib import Path
from XFELDetAna import xfelpyanatools as xana from XFELDetAna import xfelpyanatools as xana
from XFELDetAna import xfelpycaltools as xcal from XFELDetAna import xfelpycaltools as xcal
from cal_tools import h5_copy_except
from cal_tools.epix100 import epix100lib from cal_tools.epix100 import epix100lib
from cal_tools.files import DataFile
from cal_tools.tools import ( from cal_tools.tools import (
calcat_creation_time, calcat_creation_time,
get_dir_creation_date, get_dir_creation_date,
get_constant_from_db, get_constant_from_db,
load_specified_constants, load_specified_constants,
CalibrationMetadata, CalibrationMetadata,
) )
from cal_tools.step_timing import StepTimer from cal_tools.step_timing import StepTimer
from iCalibrationDB import ( from iCalibrationDB import (
Conditions, Conditions,
Constants, Constants,
) )
warnings.filterwarnings('ignore') warnings.filterwarnings('ignore')
prettyPlotting = True prettyPlotting = True
%matplotlib inline %matplotlib inline
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
x = 708 # rows of the ePix100 x = 708 # rows of the ePix100
y = 768 # columns of the ePix100 y = 768 # columns of the ePix100
if absolute_gain: if absolute_gain:
relative_gain = True relative_gain = True
plot_unit = 'ADU' plot_unit = 'ADU'
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
in_folder = Path(in_folder) in_folder = Path(in_folder)
out_folder = Path(out_folder) out_folder = Path(out_folder)
out_folder.mkdir(parents=True, exist_ok=True) out_folder.mkdir(parents=True, exist_ok=True)
run_folder = in_folder / f"r{run:04d}" run_folder = in_folder / f"r{run:04d}"
instrument_src = instrument_source_template.format( instrument_src = instrument_source_template.format(
karabo_id, receiver_template) karabo_id, receiver_template)
print(f"Correcting run: {run_folder}") print(f"Correcting run: {run_folder}")
print(f"Instrument H5File source: {instrument_src}") print(f"Instrument H5File source: {instrument_src}")
print(f"Data corrected files are stored at: {out_folder}") print(f"Data corrected files are stored at: {out_folder}")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
creation_time = calcat_creation_time(in_folder, run, creation_time) creation_time = calcat_creation_time(in_folder, run, creation_time)
print(f"Using {creation_time.isoformat()} as creation time") print(f"Using {creation_time.isoformat()} as creation time")
metadata = CalibrationMetadata(metadata_folder or out_folder) metadata = CalibrationMetadata(metadata_folder or out_folder)
# Constant paths are saved under retrieved-constants in calibration_metadata.yml. # Constant paths are saved under retrieved-constants in calibration_metadata.yml.
# NOTE: this notebook shouldn't overwrite calibration metadata file. # NOTE: this notebook shouldn't overwrite calibration metadata file.
const_yaml = metadata.get("retrieved-constants", {}) const_yaml = metadata.get("retrieved-constants", {})
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
run_dc = RunDirectory(run_folder, _use_voview=False) run_dc = RunDirectory(run_folder, _use_voview=False)
seq_files = [Path(f.filename) for f in run_dc.select(f"*{karabo_id}*").files] seq_files = [Path(f.filename) for f in run_dc.select(f"*{karabo_id}*").files]
# If a set of sequences requested to correct, # If a set of sequences requested to correct,
# adapt seq_files list. # adapt seq_files list.
if sequences != [-1]: if sequences != [-1]:
seq_files = [f for f in seq_files if any(f.match(f"*-S{s:05d}.h5") for s in sequences)] seq_files = [f for f in seq_files if any(f.match(f"*-S{s:05d}.h5") for s in sequences)]
if not len(seq_files): if not len(seq_files):
raise IndexError("No sequence files available for the selected sequences.") raise IndexError("No sequence files available for the selected sequences.")
print(f"Processing a total of {len(seq_files)} sequence files") print(f"Processing a total of {len(seq_files)} sequence files")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
step_timer = StepTimer() step_timer = StepTimer()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
step_timer.start() step_timer.start()
sensorSize = [x, y] sensorSize = [x, y]
# Sensor area will be analysed according to blocksize # Sensor area will be analysed according to blocksize
blockSize = [sensorSize[0]//2, sensorSize[1]//2] blockSize = [sensorSize[0]//2, sensorSize[1]//2]
xcal.defaultBlockSize = blockSize xcal.defaultBlockSize = blockSize
memoryCells = 1 # ePIX has no memory cells memoryCells = 1 # ePIX has no memory cells
run_parallel = False run_parallel = False
# Read control data. # Read control data.
ctrl_data = epix100lib.epix100Ctrl( ctrl_data = epix100lib.epix100Ctrl(
run_dc=run_dc, run_dc=run_dc,
instrument_src=instrument_src, instrument_src=f"{karabo_id}/DET/{receiver_template}:daqOutput",
ctrl_src=f"{karabo_id}/DET/CONTROL", ctrl_src=f"{karabo_id}/DET/CONTROL",
) )
if integration_time < 0: if integration_time < 0:
integration_time = ctrl_data.get_integration_time() integration_time = ctrl_data.get_integration_time()
integration_time_str_add = "" integration_time_str_add = ""
else: else:
integration_time_str_add = "(manual input)" integration_time_str_add = "(manual input)"
if fix_temperature < 0: if fix_temperature < 0:
temperature = ctrl_data.get_temprature() temperature = ctrl_data.get_temprature()
temperature_k = temperature + 273.15 temperature_k = temperature + 273.15
temp_str_add = "" temp_str_add = ""
else: else:
temperature_k = fix_temperature temperature_k = fix_temperature
temperature = fix_temperature - 273.15 temperature = fix_temperature - 273.15
temp_str_add = "(manual input)" temp_str_add = "(manual input)"
print(f"Bias voltage is {bias_voltage} V") print(f"Bias voltage is {bias_voltage} V")
print(f"Detector integration time is set to {integration_time} \u03BCs {integration_time_str_add}") print(f"Detector integration time is set to {integration_time} \u03BCs {integration_time_str_add}")
print(f"Mean temperature: {temperature:0.2f}°C / {temperature_k:0.2f} K {temp_str_add}") print(f"Mean temperature: {temperature:0.2f}°C / {temperature_k:0.2f} K {temp_str_add}")
print(f"Operated in vacuum: {in_vacuum}") print(f"Operated in vacuum: {in_vacuum}")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Table of sequence files to process # Table of sequence files to process
table = [(k, f) for k, f in enumerate(seq_files)] table = [(k, f) for k, f in enumerate(seq_files)]
if len(table): if len(table):
md = display(Latex(tabulate.tabulate( md = display(Latex(tabulate.tabulate(
table, table,
tablefmt='latex', tablefmt='latex',
headers=["#", "file"] headers=["#", "file"]
))) )))
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Retrieving calibration constants ## Retrieving calibration constants
As a first step, dark maps have to be loaded. As a first step, dark maps have to be loaded.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
cond_dict = { cond_dict = {
"bias_voltage": bias_voltage, "bias_voltage": bias_voltage,
"integration_time": integration_time, "integration_time": integration_time,
"temperature": temperature_k, "temperature": temperature_k,
"in_vacuum": in_vacuum, "in_vacuum": in_vacuum,
} }
dark_condition = Conditions.Dark.ePix100(**cond_dict) dark_condition = Conditions.Dark.ePix100(**cond_dict)
# update conditions with illuminated conditins. # update conditions with illuminated conditins.
cond_dict.update({ cond_dict.update({
"photon_energy": gain_photon_energy "photon_energy": gain_photon_energy
}) })
illum_condition = Conditions.Illuminated.ePix100(**cond_dict) illum_condition = Conditions.Illuminated.ePix100(**cond_dict)
const_cond = { const_cond = {
"Offset": dark_condition, "Offset": dark_condition,
"Noise": dark_condition, "Noise": dark_condition,
"RelativeGain": illum_condition, "RelativeGain": illum_condition,
} }
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
empty_constant = np.zeros((708, 768, 1), dtype=np.float32) empty_constant = np.zeros((708, 768, 1), dtype=np.float32)
if const_yaml: # Used while reproducing corrected data. if const_yaml: # Used while reproducing corrected data.
print(f"Using stored constants in {metadata.filename}") print(f"Using stored constants in {metadata.filename}")
const_data, _ = load_specified_constants(const_yaml[karabo_da]["constants"]) const_data, _ = load_specified_constants(const_yaml[karabo_da]["constants"])
for cname, cval in const_data.items(): for cname, cval in const_data.items():
if cval is None and cname != "RelativeGain": if cval is None and cname != "RelativeGain":
const_data[cname] = empty_constant const_data[cname] = empty_constant
else: # First correction attempt. else: # First correction attempt.
const_data = dict() const_data = dict()
for cname, condition in const_cond.items(): for cname, condition in const_cond.items():
# Avoid retrieving RelativeGain, if not needed for correction. # Avoid retrieving RelativeGain, if not needed for correction.
if cname == "RelativeGain" and not relative_gain: if cname == "RelativeGain" and not relative_gain:
const_data[cname] = None const_data[cname] = None
else: else:
const_data[cname] = get_constant_from_db( const_data[cname] = get_constant_from_db(
karabo_id=karabo_id, karabo_id=karabo_id,
karabo_da=karabo_da, karabo_da=karabo_da,
constant=getattr(Constants.ePix100, cname)(), constant=getattr(Constants.ePix100, cname)(),
condition=condition, condition=condition,
empty_constant=None if cname == "RelativeGain" else empty_constant, empty_constant=None if cname == "RelativeGain" else empty_constant,
cal_db_interface=cal_db_interface, cal_db_interface=cal_db_interface,
creation_time=creation_time, creation_time=creation_time,
print_once=2, print_once=2,
timeout=cal_db_timeout timeout=cal_db_timeout
) )
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if relative_gain and const_data.get("RelativeGain", None) is None: if relative_gain and const_data.get("RelativeGain", None) is None:
print( print(
"WARNING: RelativeGain map is requested, but not found.\n" "WARNING: RelativeGain map is requested, but not found.\n"
"No gain correction will be applied" "No gain correction will be applied"
) )
relative_gain = False relative_gain = False
absolute_gain = False absolute_gain = False
# Initializing some parameters. # Initializing some parameters.
hscale = 1 hscale = 1
stats = True stats = True
hrange = np.array([-50, 1000]) hrange = np.array([-50, 1000])
nbins = hrange[1] - hrange[0] nbins = hrange[1] - hrange[0]
commonModeBlockSize = [x//2, y//2] commonModeBlockSize = [x//2, y//2]
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
histCalOffsetCor = xcal.HistogramCalculator( histCalOffsetCor = xcal.HistogramCalculator(
sensorSize, sensorSize,
bins=nbins, bins=nbins,
range=hrange, range=hrange,
parallel=run_parallel, parallel=run_parallel,
nCells=memoryCells, nCells=memoryCells,
blockSize=blockSize blockSize=blockSize
) )
# *****************Histogram Calculators****************** # # *****************Histogram Calculators****************** #
histCalCor = xcal.HistogramCalculator( histCalCor = xcal.HistogramCalculator(
sensorSize, sensorSize,
bins=1050, bins=1050,
range=[-50, 1000], range=[-50, 1000],
parallel=run_parallel, parallel=run_parallel,
nCells=memoryCells, nCells=memoryCells,
blockSize=blockSize blockSize=blockSize
) )
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if common_mode: if common_mode:
histCalCMCor = xcal.HistogramCalculator( histCalCMCor = xcal.HistogramCalculator(
sensorSize, sensorSize,
bins=nbins, bins=nbins,
range=hrange, range=hrange,
parallel=run_parallel, parallel=run_parallel,
nCells=memoryCells, nCells=memoryCells,
blockSize=blockSize, blockSize=blockSize,
) )
cmCorrectionB = xcal.CommonModeCorrection( cmCorrectionB = xcal.CommonModeCorrection(
shape=sensorSize, shape=sensorSize,
blockSize=commonModeBlockSize, blockSize=commonModeBlockSize,
orientation='block', orientation='block',
nCells=memoryCells, nCells=memoryCells,
noiseMap=const_data['Noise'], noiseMap=const_data['Noise'],
runParallel=run_parallel, runParallel=run_parallel,
parallel=run_parallel, parallel=run_parallel,
stats=stats, stats=stats,
minFrac=cm_min_frac, minFrac=cm_min_frac,
noiseSigma=cm_noise_sigma, noiseSigma=cm_noise_sigma,
) )
cmCorrectionR = xcal.CommonModeCorrection( cmCorrectionR = xcal.CommonModeCorrection(
shape=sensorSize, shape=sensorSize,
blockSize=commonModeBlockSize, blockSize=commonModeBlockSize,
orientation='row', orientation='row',
nCells=memoryCells, nCells=memoryCells,
noiseMap=const_data['Noise'], noiseMap=const_data['Noise'],
runParallel=run_parallel, runParallel=run_parallel,
parallel=run_parallel, parallel=run_parallel,
stats=stats, stats=stats,
minFrac=cm_min_frac, minFrac=cm_min_frac,
noiseSigma=cm_noise_sigma, noiseSigma=cm_noise_sigma,
) )
cmCorrectionC = xcal.CommonModeCorrection( cmCorrectionC = xcal.CommonModeCorrection(
shape=sensorSize, shape=sensorSize,
blockSize=commonModeBlockSize, blockSize=commonModeBlockSize,
orientation='col', orientation='col',
nCells=memoryCells, nCells=memoryCells,
noiseMap=const_data['Noise'], noiseMap=const_data['Noise'],
runParallel=run_parallel, runParallel=run_parallel,
parallel=run_parallel, parallel=run_parallel,
stats=stats, stats=stats,
minFrac=cm_min_frac, minFrac=cm_min_frac,
noiseSigma=cm_noise_sigma, noiseSigma=cm_noise_sigma,
) )
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if relative_gain: if relative_gain:
gain_cnst = np.median(const_data["RelativeGain"]) gain_cnst = np.median(const_data["RelativeGain"])
hscale = gain_cnst hscale = gain_cnst
plot_unit = 'keV' plot_unit = 'keV'
if photon_energy > 0: if photon_energy > 0:
plot_unit = '$\gamma$' plot_unit = '$\gamma$'
hscale /= photon_energy hscale /= photon_energy
gainCorrection = xcal.RelativeGainCorrection( gainCorrection = xcal.RelativeGainCorrection(
sensorSize, sensorSize,
gain_cnst/const_data["RelativeGain"][..., None], gain_cnst/const_data["RelativeGain"][..., None],
nCells=memoryCells, nCells=memoryCells,
parallel=run_parallel, parallel=run_parallel,
blockSize=blockSize, blockSize=blockSize,
gains=None, gains=None,
) )
histCalRelGainCor = xcal.HistogramCalculator( histCalRelGainCor = xcal.HistogramCalculator(
sensorSize, sensorSize,
bins=nbins, bins=nbins,
range=hrange, range=hrange,
parallel=run_parallel, parallel=run_parallel,
nCells=memoryCells, nCells=memoryCells,
blockSize=blockSize blockSize=blockSize
) )
if absolute_gain: if absolute_gain:
histCalAbsGainCor = xcal.HistogramCalculator( histCalAbsGainCor = xcal.HistogramCalculator(
sensorSize, sensorSize,
bins=nbins, bins=nbins,
range=hrange*hscale, range=hrange*hscale,
parallel=run_parallel, parallel=run_parallel,
nCells=memoryCells, nCells=memoryCells,
blockSize=blockSize blockSize=blockSize
) )
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if pattern_classification : if pattern_classification :
patternClassifier = xcal.PatternClassifier( patternClassifier = xcal.PatternClassifier(
[x, y], [x, y],
const_data["Noise"], const_data["Noise"],
split_evt_primary_threshold, split_evt_primary_threshold,
split_evt_secondary_threshold, split_evt_secondary_threshold,
split_evt_mip_threshold, split_evt_mip_threshold,
tagFirstSingles=0, tagFirstSingles=0,
nCells=memoryCells, nCells=memoryCells,
allowElongated=False, allowElongated=False,
blockSize=[x, y], blockSize=[x, y],
parallel=run_parallel, parallel=run_parallel,
) )
histCalCSCor = xcal.HistogramCalculator( histCalCSCor = xcal.HistogramCalculator(
sensorSize, sensorSize,
bins=nbins, bins=nbins,
range=hrange, range=hrange,
parallel=run_parallel, parallel=run_parallel,
nCells=memoryCells, nCells=memoryCells,
blockSize=blockSize, blockSize=blockSize,
) )
histCalGainCorClusters = xcal.HistogramCalculator( histCalGainCorClusters = xcal.HistogramCalculator(
sensorSize, sensorSize,
bins=nbins, bins=nbins,
range=hrange*hscale, range=hrange*hscale,
parallel=run_parallel, parallel=run_parallel,
nCells=memoryCells, nCells=memoryCells,
blockSize=blockSize blockSize=blockSize
) )
histCalGainCorSingles = xcal.HistogramCalculator( histCalGainCorSingles = xcal.HistogramCalculator(
sensorSize, sensorSize,
bins=nbins, bins=nbins,
range=hrange*hscale, range=hrange*hscale,
parallel=run_parallel, parallel=run_parallel,
nCells=memoryCells, nCells=memoryCells,
blockSize=blockSize blockSize=blockSize
) )
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Applying corrections ## Applying corrections
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def correct_train(wid, index, tid, d): def correct_train(wid, index, tid, d):
d = d[..., np.newaxis].astype(np.float32) d = d[pixel_data[0]][pixel_data[1]][..., np.newaxis].astype(np.float32)
d = np.compress( d = np.compress(
np.any(d > 0, axis=(0, 1)), d, axis=2) np.any(d > 0, axis=(0, 1)), d, axis=2)
# Offset correction. # Offset correction.
d -= const_data["Offset"] d -= const_data["Offset"]
histCalOffsetCor.fill(d) histCalOffsetCor.fill(d)
# Common Mode correction. # Common Mode correction.
if common_mode: if common_mode:
# Block CM # Block CM
d = cmCorrectionB.correct(d) d = cmCorrectionB.correct(d)
# Row CM # Row CM
d = cmCorrectionR.correct(d) d = cmCorrectionR.correct(d)
# COL CM # COL CM
d = cmCorrectionC.correct(d) d = cmCorrectionC.correct(d)
histCalCMCor.fill(d) histCalCMCor.fill(d)
# relative gain correction. # relative gain correction.
if relative_gain: if relative_gain:
d = gainCorrection.correct(d) d = gainCorrection.correct(d)
histCalRelGainCor.fill(d) histCalRelGainCor.fill(d)
"""The gain correction is currently applying """The gain correction is currently applying
an absolute correction (not a relative correction an absolute correction (not a relative correction
as the implied by the name); as the implied by the name);
it changes the scale (the unit of measurement) it changes the scale (the unit of measurement)
of the data from ADU to either keV or n_of_photons. of the data from ADU to either keV or n_of_photons.
But the pattern classification relies on comparing But the pattern classification relies on comparing
data with the noise map, which is still in ADU. data with the noise map, which is still in ADU.
The best solution is to do a relative gain The best solution is to do a relative gain
correction first and apply the global absolute correction first and apply the global absolute
gain to the data at the end, after clustering. gain to the data at the end, after clustering.
""" """
if pattern_classification: if pattern_classification:
d_clu, patterns = patternClassifier.classify(d) d_clu, patterns = patternClassifier.classify(d)
d_clu[d_clu < (split_evt_primary_threshold*const_data["Noise"])] = 0 d_clu[d_clu < (split_evt_primary_threshold*const_data["Noise"])] = 0
data_clu[index, ...] = np.squeeze(d_clu) data_clu[index, ...] = np.squeeze(d_clu)
data_patterns[index, ...] = np.squeeze(patterns) data_patterns[index, ...] = np.squeeze(patterns)
histCalCSCor.fill(d_clu) histCalCSCor.fill(d_clu)
# absolute gain correction # absolute gain correction
# changes data from ADU to keV (or n. of photons) # changes data from ADU to keV (or n. of photons)
if absolute_gain: if absolute_gain:
d = d * gain_cnst d = d * gain_cnst
if photon_energy > 0: if photon_energy > 0:
d /= photon_energy d /= photon_energy
histCalAbsGainCor.fill(d) histCalAbsGainCor.fill(d)
if pattern_classification: if pattern_classification:
# Modify pattern classification. # Modify pattern classification.
d_clu = d_clu * gain_cnst d_clu = d_clu * gain_cnst
if photon_energy > 0: if photon_energy > 0:
d_clu /= photon_energy d_clu /= photon_energy
data_clu[index, ...] = np.squeeze(d_clu) data_clu[index, ...] = np.squeeze(d_clu)
histCalGainCorClusters.fill(d_clu) histCalGainCorClusters.fill(d_clu)
d_sing = d_clu[patterns==100] # pattern 100 corresponds to single photons events d_sing = d_clu[patterns==100] # pattern 100 corresponds to single photons events
if len(d_sing): if len(d_sing):
histCalGainCorSingles.fill(d_sing) histCalGainCorSingles.fill(d_sing)
data[index, ...] = np.squeeze(d) data[index, ...] = np.squeeze(d)
histCalCor.fill(d) histCalCor.fill(d)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
pixel_data = (instrument_src, "data.image.pixels")
# 10 is a number chosen after testing 1 ... 71 parallel threads # 10 is a number chosen after testing 1 ... 71 parallel threads
context = psh.context.ThreadContext(num_workers=10) context = psh.context.ThreadContext(num_workers=10)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
empty_seq = 0
for f in seq_files: for f in seq_files:
seq_dc = H5File(f) seq_dc = H5File(f)
# Save corrected data in an output file with name
# of corresponding raw sequence file. n_imgs = seq_dc.get_data_counts(*pixel_data).shape[0]
out_file = out_folder / f.name.replace("RAW", "CORR")
# Data shape in seq_dc excluding trains with empty images. # Data shape in seq_dc excluding trains with empty images.
ishape = seq_dc[instrument_src, "data.image.pixels"].shape dshape = seq_dc[pixel_data].shape
corr_ntrains = ishape[0] dataset_chunk = ((chunk_size_idim,) + dshape[1:]) # e.g. (1, pixels_x, pixels_y)
all_train_ids = seq_dc.train_ids
if n_imgs - dshape[0] != 0:
# Raise a WARNING if this sequence has no trains to correct. print(f"- WARNING: {f} has {n_imgs - dshape[0]} trains with empty data.")
# Otherwise, print number of trains with no data.
if corr_ntrains == 0:
warning(f"No trains to correct for {f.name}: "
"Skipping the processing of this file.")
empty_seq += 1
continue
elif len(all_train_ids) != corr_ntrains:
print(f"{f.name} has {len(all_train_ids) - corr_ntrains} trains with missing data.")
# This parameter is only used for testing. # This parameter is only used for testing.
if limit_trains > 0: if limit_images > 0:
print(f"\nCorrected trains are limited to: {limit_trains} trains") n_imgs = min(n_imgs, limit_images)
corr_ntrains = min(corr_ntrains, limit_trains)
oshape = (corr_ntrains, *ishape[1:])
data = context.alloc(shape=oshape, dtype=np.float32) data = context.alloc(shape=dshape, dtype=np.float32)
if pattern_classification: if pattern_classification:
data_clu = context.alloc(shape=oshape, dtype=np.float32) data_clu = context.alloc(shape=dshape, dtype=np.float32)
data_patterns = context.alloc(shape=oshape, dtype=np.int32) data_patterns = context.alloc(shape=dshape, dtype=np.int32)
step_timer.start() # Correct data.
# Overwrite seq_dc after eliminating empty trains or/and applying limited images.
seq_dc = seq_dc.select(
instrument_src, "*", require_all=True).select_trains(np.s_[:corr_ntrains])
pixel_data = seq_dc[instrument_src, "data.image.pixels"]
context.map(correct_train, pixel_data)
step_timer.done_step(f'Correcting {corr_ntrains} trains.')
step_timer.start()
step_timer.start() # Write corrected data. context.map(
correct_train, seq_dc.select(
# Create CORR files and add corrected data sections. *pixel_data, require_all=True).select_trains(np.s_[:n_imgs])
image_counts = seq_dc[instrument_src, "data.image.pixels"].data_counts(labelled=False) )
step_timer.done_step(f'Correcting {n_imgs} trains.')
# Write corrected data.
with DataFile(out_file, 'w') as ofile:
dataset_chunk = ((chunk_size_idim,) + oshape[1:]) # e.g. (1, pixels_x, pixels_y)
# Create INDEX datasets.
ofile.create_index(seq_dc.train_ids, from_file=seq_dc.files[0])
# Create Instrument section to later add corrected datasets.
outp_source = ofile.create_instrument_source(instrument_src)
# Create count/first datasets at INDEX source.
outp_source.create_index(data=image_counts)
# Store uncorrected RAW image datasets for the corrected trains.
data_raw_fields = [ # /data/
'ambTemp', 'analogCurr', 'analogInputVolt', 'backTemp',
'digitalInputVolt', 'guardCurr', 'relHumidity',
]
for field in data_raw_fields:
field_arr = seq_dc[instrument_src, f"data.{field}"].ndarray()
outp_source.create_key(
f"data.{field}", data=field_arr,
chunks=(chunk_size_idim, *field_arr.shape[1:]))
image_raw_fields = [ # /data/image/ # Store detector h5 information in the corrected file
'binning', 'bitsPerPixel', 'dimTypes', 'dims', # and deselect data to correct and store later.
'encoding', 'flipX', 'flipY', 'roiOffsets', 'rotation', step_timer.start()
]
for field in image_raw_fields:
field_arr = seq_dc[instrument_src, f"data.image.{field}"].ndarray()
outp_source.create_key( out_file = out_folder / f.name.replace("RAW", "CORR")
f"data.image.{field}", data=field_arr, data_path = "INSTRUMENT/"+instrument_src+"/data/image"
chunks=(chunk_size_idim, *field_arr.shape[1:])) pixels_path = f"{data_path}/pixels"
# Add main corrected `data.image.pixels` dataset and store corrected data. # First copy all raw data source to the corrected file,
outp_source.create_key( # while excluding the raw data image /data/image/pixels.
"data.image.pixels", data=data, chunks=dataset_chunk) with h5py.File(out_file, 'w') as ofile:
# Copy RAW non-calibrated sources.
with h5py.File(f, 'r') as sfile:
h5_copy_except.h5_copy_except_paths(
sfile, ofile,
[pixels_path])
# Create dataset in CORR h5 file and add corrected images.
dataset = ofile.create_dataset(
pixels_path,
data=data,
chunks=dataset_chunk,
dtype=np.float32)
if pattern_classification: if pattern_classification:
# Add main corrected `data.image.pixels` dataset and store corrected data. # Save /data/image/pixels_classified in corrected file.
outp_source.create_key( datasetc = ofile.create_dataset(
"data.image.pixels_classified", data=data_clu, chunks=dataset_chunk) f"{data_path}/pixels_classified",
outp_source.create_key( data=data_clu,
"data.image.patterns", data=data_clu, chunks=dataset_chunk) chunks=dataset_chunk,
dtype=np.float32)
# Create METDATA datasets
ofile.create_metadata(like=seq_dc) # Save /data/image/patterns in corrected file.
datasetp = ofile.create_dataset(
f"{data_path}/patterns",
data=data_patterns,
chunks=dataset_chunk,
dtype=np.int32)
step_timer.done_step('Storing data.') step_timer.done_step('Storing data.')
if empty_seq == len(seq_files):
warning("No valid trains for RAW data to correct.")
exit(0)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
ho, eo, co, so = histCalCor.get() ho, eo, co, so = histCalCor.get()
d = [{ d = [{
'x': co, 'x': co,
'y': ho, 'y': ho,
'y_err': np.sqrt(ho[:]), 'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid', 'drawstyle': 'steps-mid',
'errorstyle': 'bars', 'errorstyle': 'bars',
'errorcoarsing': 2, 'errorcoarsing': 2,
'label': 'Total corr.' 'label': 'Total corr.'
}] }]
ho, eo, co, so = histCalOffsetCor.get() ho, eo, co, so = histCalOffsetCor.get()
d.append({ d.append({
'x': co, 'x': co,
'y': ho, 'y': ho,
'y_err': np.sqrt(ho[:]), 'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid', 'drawstyle': 'steps-mid',
'errorstyle': 'bars', 'errorstyle': 'bars',
'errorcoarsing': 2, 'errorcoarsing': 2,
'label': 'Offset corr.' 'label': 'Offset corr.'
}) })
if common_mode: if common_mode:
ho, eo, co, so = histCalCMCor.get() ho, eo, co, so = histCalCMCor.get()
d.append({ d.append({
'x': co, 'x': co,
'y': ho, 'y': ho,
'y_err': np.sqrt(ho[:]), 'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid', 'drawstyle': 'steps-mid',
'errorstyle': 'bars', 'errorstyle': 'bars',
'errorcoarsing': 2, 'errorcoarsing': 2,
'label': 'CM corr.' 'label': 'CM corr.'
}) })
if relative_gain : if relative_gain :
ho, eo, co, so = histCalRelGainCor.get() ho, eo, co, so = histCalRelGainCor.get()
d.append({ d.append({
'x': co, 'x': co,
'y': ho, 'y': ho,
'y_err': np.sqrt(ho[:]), 'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid', 'drawstyle': 'steps-mid',
'errorstyle': 'bars', 'errorstyle': 'bars',
'errorcoarsing': 2, 'errorcoarsing': 2,
'label': 'Relative gain corr.' 'label': 'Relative gain corr.'
}) })
if pattern_classification: if pattern_classification:
ho, eo, co, so = histCalCSCor.get() ho, eo, co, so = histCalCSCor.get()
d.append({ d.append({
'x': co, 'x': co,
'y': ho, 'y': ho,
'y_err': np.sqrt(ho[:]), 'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid', 'drawstyle': 'steps-mid',
'errorstyle': 'bars', 'errorstyle': 'bars',
'errorcoarsing': 2, 'errorcoarsing': 2,
'label': 'Charge sharing corr.' 'label': 'Charge sharing corr.'
}) })
fig = xana.simplePlot( fig = xana.simplePlot(
d, aspect=1, x_label=f'Energy (ADU)', d, aspect=1, x_label=f'Energy (ADU)',
y_label='Number of occurrences', figsize='2col', y_label='Number of occurrences', figsize='2col',
y_log=True, x_range=(-50, 500), y_log=True, x_range=(-50, 500),
legend='top-center-frame-2col', legend='top-center-frame-2col',
) )
plt.title(f'run {run} - {karabo_da}') plt.title(f'run {run} - {karabo_da}')
plt.grid() plt.grid()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if absolute_gain : if absolute_gain :
d=[] d=[]
ho, eo, co, so = histCalAbsGainCor.get() ho, eo, co, so = histCalAbsGainCor.get()
d.append({ d.append({
'x': co, 'x': co,
'y': ho, 'y': ho,
'y_err': np.sqrt(ho[:]), 'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid', 'drawstyle': 'steps-mid',
'errorstyle': 'bars', 'errorstyle': 'bars',
'errorcoarsing': 2, 'errorcoarsing': 2,
'label': 'Absolute gain corr.' 'label': 'Absolute gain corr.'
}) })
if pattern_classification: if pattern_classification:
ho, eo, co, so = histCalGainCorClusters.get() ho, eo, co, so = histCalGainCorClusters.get()
d.append({ d.append({
'x': co, 'x': co,
'y': ho, 'y': ho,
'y_err': np.sqrt(ho[:]), 'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid', 'drawstyle': 'steps-mid',
'errorstyle': 'bars', 'errorstyle': 'bars',
'errorcoarsing': 2, 'errorcoarsing': 2,
'label': 'Charge sharing corr.' 'label': 'Charge sharing corr.'
}) })
ho, eo, co, so = histCalGainCorSingles.get() ho, eo, co, so = histCalGainCorSingles.get()
d.append({ d.append({
'x': co, 'x': co,
'y': ho, 'y': ho,
'y_err': np.sqrt(ho[:]), 'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid', 'drawstyle': 'steps-mid',
'errorstyle': 'bars', 'errorstyle': 'bars',
'errorcoarsing': 2, 'errorcoarsing': 2,
'label': 'Isolated photons (singles)' 'label': 'Isolated photons (singles)'
}) })
fig = xana.simplePlot( fig = xana.simplePlot(
d, aspect=1, x_label=f'Energy ({plot_unit})', d, aspect=1, x_label=f'Energy ({plot_unit})',
y_label='Number of occurrences', figsize='2col', y_label='Number of occurrences', figsize='2col',
y_log=True, y_log=True,
x_range=np.array((-50, 500))*hscale, x_range=np.array((-50, 500))*hscale,
legend='top-center-frame-2col', legend='top-center-frame-2col',
) )
plt.grid() plt.grid()
plt.title(f'run {run} - {karabo_da}') plt.title(f'run {run} - {karabo_da}')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Mean Image of the corrected data ## Mean Image of the corrected data
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
step_timer.start() step_timer.start()
fig = xana.heatmapPlot( fig = xana.heatmapPlot(
np.nanmedian(data, axis=0), np.nanmedian(data, axis=0),
x_label='Columns', y_label='Rows', x_label='Columns', y_label='Rows',
lut_label=f'Signal ({plot_unit})', lut_label=f'Signal ({plot_unit})',
x_range=(0, y), x_range=(0, y),
y_range=(0, x), y_range=(0, x),
vmin=-50, vmax=50) vmin=-50, vmax=50)
step_timer.done_step(f'Plotting mean image of {data.shape[0]} trains.') step_timer.done_step(f'Plotting mean image of {data.shape[0]} trains.')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Single Shot of the corrected data ## Single Shot of the corrected data
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
step_timer.start() step_timer.start()
fig = xana.heatmapPlot( fig = xana.heatmapPlot(
data[0, ...], data[0, ...],
x_label='Columns', y_label='Rows', x_label='Columns', y_label='Rows',
lut_label=f'Signal ({plot_unit})', lut_label=f'Signal ({plot_unit})',
x_range=(0, y), x_range=(0, y),
y_range=(0, x), y_range=(0, x),
vmin=-50, vmax=50) vmin=-50, vmax=50)
step_timer.done_step(f'Plotting single shot of corrected data.') step_timer.done_step(f'Plotting single shot of corrected data.')
``` ```
......
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