Skip to content
Snippets Groups Projects

[PNCCD][CORRECT] Using calcat interface

Merged [PNCCD][CORRECT] Using calcat interface
feat/pnccd_calcat_interface into master
All threads resolved!
Merged Karim Ahmed requested to merge feat/pnccd_calcat_interface into master
All threads resolved!
Viewing commit a30fac59
Show latest version
2 files
+ 129
− 113
Compare changes
  • Side-by-side
  • Inline
Files
2
notebooks/pnCCD/Correct_pnCCD_NBC.ipynb
+ 65
− 62
%% Cell type:markdown id: tags:
# pnCCD Data Correction #
Authors: DET Group, Modified by Kiana Setoodehnia - Version 5.0
The following notebook provides offset, common mode, relative gain, split events and pattern classification corrections of images acquired with the pnCCD. This notebook *does not* yet correct for charge transfer inefficiency.
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/SQS/202031/p900166/raw" # input folder
out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/remove/pnccd_correct" # output folder
metadata_folder = "" # Directory containing calibration_metadata.yml when run by xfel-calibrate
run = 347 # which run to read data from
sequences = [-1] # sequences to correct, set to -1 for all, range allowed
sequences_per_node = 1 # number of sequences running on the same slurm node.
karabo_da = 'PNCCD01' # data aggregators
karabo_id = "SQS_NQS_PNCCD1MP" # karabo prefix of PNCCD devices
receiver_id = "PNCCD_FMT-0" # inset for receiver devices
path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data
instrument_source_template = '{}/CAL/{}:output' # template for data source name, will be filled with karabo_id and receiver_id.
# Parameters affecting data correction.
commonModeAxis = 0 # axis along which common mode will be calculated, 0 = row, and 1 = column
commonModeBlockSize = [512, 512] # size of the detector in pixels for common mode calculations
split_evt_primary_threshold = 4. # primary threshold for split event classification in terms of n sigma noise
split_evt_secondary_threshold = 3. # secondary threshold for split event classification in terms of n sigma noise
saturated_threshold = 32000. # full well capacity in ADU
# Conditions for retrieving calibration constants
fix_temperature_top = 0. # fix temperature for top sensor in K, set to 0. to use value from slow data.
fix_temperature_bot = 0. # fix temperature for bottom senspr in K, set to 0. to use value from slow data.
gain = -1 # the detector's gain setting. Set to -1 to use the value from the slow data.
bias_voltage = 0. # the detector's bias voltage. set to 0. to use value from slow data.
integration_time = 70 # detector's integration time
photon_energy = 1.6 # Al fluorescence in keV
# Parameters for the calibration database.
cal_db_interface = "tcp://max-exfl016:8015" # calibration DB interface to use
cal_db_timeout = 300000 # timeout on caldb requests
creation_time = "" # The timestamp to use with Calibration DB. Required Format: "YYYY-MM-DD hh:mm:ss" e.g. 2019-07-04 11:02:41
# Booleans for selecting corrections to apply.
only_offset = False # Only, apply offset.
common_mode = True # Apply common mode correction
relgain = True # Apply relative gain correction
pattern_classification = True # classify split events
# parameters affecting stored output data.
chunk_size_idim = 1 # H5 chunking size of output data
# ONLY FOR TESTING
limit_images = 0 # this parameter is used for limiting number of images to correct from a sequence file. ONLY FOR TESTING.
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
# Here the herarichy and dependability for correction booleans are defined
corr_bools = {}
corr_bools["only_offset"] = only_offset
# Apply offset only.
if not only_offset:
corr_bools["relgain"] = relgain
corr_bools["common_mode"] = common_mode
corr_bools["pattern_class"] = pattern_classification
```
%% Cell type:code id: tags:
``` python
import datetime
import os
import warnings
from logging import warning
from pathlib import Path
warnings.filterwarnings('ignore')
import h5py
import matplotlib.pyplot as plt
import numpy as np
import pasha as psh
from IPython.display import Markdown, display
from extra_data import H5File, RunDirectory
from prettytable import PrettyTable
%matplotlib inline
from calibration_client import CalibrationClient
from XFELDetAna import xfelpyanatools as xana
from XFELDetAna import xfelpycaltools as xcal
from cal_tools import pnccdlib
from cal_tools.calcat_interface import PNCCD_CalibrationData
from cal_tools.restful_config import restful_config
from cal_tools.tools import (
calcat_creation_time,
get_dir_creation_date,
get_constant_from_db_and_time,
get_random_db_interface,
load_specified_constants,
load_constants_dict,
CalibrationMetadata,
)
from cal_tools.step_timing import StepTimer
from cal_tools import h5_copy_except
from iCalibrationDB import Conditions, Constants
from iCalibrationDB.detectors import DetectorTypes
```
%% Cell type:code id: tags:
``` python
# Calibration Database Settings, and Some Initial Run Parameters & Paths:
display(Markdown('### Initial Settings and Paths'))
# Sensor size and block size definitions (important for common mode and other corrections):
pixels_x = 1024 # rows of pnCCD in pixels
pixels_y = 1024 # columns of pnCCD in pixels
in_folder = Path(in_folder)
sensorSize = [pixels_x, pixels_y]
# For xcal.HistogramCalculators.
blockSize = [pixels_x//2, pixels_y//2] # sensor area will be analysed according to blockSize.
print(f"pnCCD size is: {pixels_x}x{pixels_y} pixels.")
print(f'Calibration database interface selected: {cal_db_interface}')
# Paths to the data:
instrument_src = instrument_source_template.format(karabo_id, receiver_id)
print(f"Instrument H5File source: {instrument_src}\n")
# Run's creation time:
creation_time = calcat_creation_time(in_folder, run, creation_time)
print(f"Creation time: {creation_time}")
```
%% Cell type:code id: tags:
``` python
step_timer = StepTimer()
```
%% Cell type:code id: tags:
``` python
run_dc = RunDirectory(in_folder / f"r{run:04d}", _use_voview=False)
# Output Folder Creation:
os.makedirs(out_folder, exist_ok=True)
# NOTE: this notebook shouldn't overwrite calibration metadata file.
metadata = CalibrationMetadata(metadata_folder or out_folder)
# Constant paths are saved under retrieved-constants in calibration_metadata.yml
const_yaml = metadata.get("retrieved-constants", {})
# extract control data
step_timer.start()
ctrl_data = pnccdlib.PnccdCtrl(run_dc, karabo_id)
if bias_voltage == 0.:
bias_voltage = ctrl_data.get_bias_voltage()
if gain == -1:
gain = ctrl_data.get_gain()
if fix_temperature_top == 0:
fix_temperature_top = ctrl_data.get_fix_temperature_top()
if fix_temperature_bot == 0:
fix_temperature_bot = ctrl_data.get_fix_temperature_bot()
step_timer.done_step("Reading control parameters.")
# Printing the Parameters Read from the Data File:
display(Markdown('### Detector Parameters'))
print(f"Bias voltage is {bias_voltage:0.1f} V.")
print(f"Detector gain is set to 1/{int(gain)}.")
print(f"Detector integration time is set to {integration_time} ms")
print(f"Top pnCCD sensor is at temperature of {fix_temperature_top:0.2f} K")
print(f"Bottom pnCCD sensor is at temperature of {fix_temperature_bot:0.2f} K")
```
%% Cell type:code id: tags:
``` python
seq_files = []
for f in run_dc.select(instrument_src).files:
fpath = Path(f.filename)
if fpath.match(f"*{karabo_da}*.h5"):
seq_files.append(fpath)
if sequences != [-1]:
seq_files = sorted([f for f in seq_files if any(f.match(f"*-S{s:05d}.h5") for s in sequences)])
print(f"Processing a total of {len(seq_files)} sequence files:")
print(*seq_files, sep='\n')
```
%% Cell type:code id: tags:
``` python
gain_k = [k for k, v in pnccdlib.VALID_GAINS.items() if v == gain][0]
if gain_k == 'a':
split_evt_mip_threshold = 1000. # MIP threshold in ADU for event classification (10 times average noise)
# Each xcal.HistogramCalculator requires a total number of bins and a binning range. We define these
# using a dictionary:
# For all xcal histograms:
Hist_Bin_Dict = {
"bins": 35000, # number of bins
"bin_range": [0, 35000]
}
# For the numpy histograms on the last cell of the notebook:
Event_Bin_Dict = {
"event_bins": 1000, # number of bins
"b_range": [0, 35000] # bin range
}
elif gain_k == 'b':
split_evt_mip_threshold = 270. # 10 times the average noise
Hist_Bin_Dict = {
"bins": 10000,
"bin_range": [0, 10000]
}
Event_Bin_Dict = {
"event_bins": 1000,
"b_range": [0, 10000]
}
elif gain_k == 'c':
split_evt_mip_threshold = 110. # 10 times the average noise
Hist_Bin_Dict = {
"bins": 3000,
"bin_range": [0, 3000]
}
Event_Bin_Dict = {
"event_bins": 1000,
"b_range": [0, 3000]
}
elif gain_k == 'd':
split_evt_mip_threshold = 90. # 10 times the average noise
Hist_Bin_Dict = {
"bins": 1000,
"bin_range": [0, 1000]
}
Event_Bin_Dict = {
"event_bins": 1000,
"b_range": [0, 1000]
}
elif gain_k == 'e':
split_evt_mip_threshold = 90. # 10 times the average noise
Hist_Bin_Dict = {
"bins": 500,
"bin_range": [0, 500]
}
Event_Bin_Dict = {
"event_bins": 500,
"b_range": [0, 500]
}
else:
split_evt_mip_threshold = 90. # 10 times the average noise
Hist_Bin_Dict = {
"bins": 220,
"bin_range": [0, 220]
}
Event_Bin_Dict = {
"event_bins": 220,
"b_range": [0, 220]
}
bins = Hist_Bin_Dict["bins"]
bin_range = Hist_Bin_Dict["bin_range"]
event_bins = Event_Bin_Dict["event_bins"]
b_range = Event_Bin_Dict["b_range"]
```
%% Cell type:code id: tags:
``` python
# Connect to CalCat.
calcat_config = restful_config['calcat']
client = CalibrationClient(
base_api_url=calcat_config['base-api-url'],
use_oauth2=calcat_config['use-oauth2'],
client_id=calcat_config['user-id'],
client_secret=calcat_config['user-secret'],
user_email=calcat_config['user-email'],
token_url=calcat_config['token-url'],
refresh_url=calcat_config['refresh-url'],
auth_url=calcat_config['auth-url'],
scope='')
```
%% Cell type:markdown id: tags:
As a first step, dark constants have to be retrieved from the calibration database
%% Cell type:code id: tags:
``` python
display(Markdown("### Constants retrieval"))
step_timer.start()
conditions_dict = {
"bias_voltage": bias_voltage,
"integration_time": integration_time,
"gain_setting": gain,
"temperature": fix_temperature_top,
"pixels_x": pixels_x,
"pixels_y": pixels_y,
}
# Dark condition
dark_condition = Conditions.Dark.CCD(**conditions_dict)
# Add photon energy.
conditions_dict.update({"photon_energy": photon_energy})
illum_condition = Conditions.Illuminated.CCD(**conditions_dict)
# A dictionary for initializing constants. {cname: empty constant array}
empty_constants = {
"Offset": np.zeros((pixels_x, pixels_y, 1), dtype=np.float32),
"Noise": np.zeros((pixels_x, pixels_y, 1), dtype=np.float32),
"BadPixelsDark": np.zeros((pixels_x, pixels_y, 1), dtype=np.uint32),
"RelativeGain": np.zeros((pixels_x, pixels_y), dtype=np.float32),
}
constant_names = ["OffsetCCD", "NoiseCCD", "BadPixelsDarkCCD"]
if relgain:
constant_names += ["RelativeGainCCD"]
if const_yaml: # Used while reproducing corrected data.
print(f"Using stored constants in {metadata.filename}")
constants, when = load_specified_constants(
const_yaml[karabo_da]["constants"], empty_constants
)
constants, _ = load_constants_dict(const_yaml[karabo_da]["constants"])
else:
constants = dict()
when = dict()
for cname, cempty in empty_constants.items():
# No need for retrieving RelativeGain, if not used for correction.
if not corr_bools.get("relgain") and cname == "RelativeGain":
continue
constants[cname], when[cname] = get_constant_from_db_and_time(
karabo_id,
karabo_da,
constant=getattr(Constants.CCD(DetectorTypes.pnCCD), cname)(),
condition=illum_condition if cname == "RelativeGain" else dark_condition,
empty_constant=cempty,
cal_db_interface=get_random_db_interface(cal_db_interface),
creation_time=creation_time,
)
pnccd_cal = PNCCD_CalibrationData(
detector_name=karabo_id,
sensor_bias_voltage=bias_voltage,
integration_time=integration_time,
sensor_temperature=fix_temperature_top,
gain_setting=gain,
event_at=creation_time,
snapshot_at=creation_time,
client=client,
)
constants = pnccd_cal.ndarray_map(calibrations=constant_names)[karabo_da]
# Validate the constants availability and raise/warn correspondingly.
missing_dark_constants = set(
c for c in ["OffsetCCD", "NoiseCCD", "BadPixelsDarkCCD"] if c not in constants.keys())
if missing_dark_constants:
raise KeyError(
f"Dark constants {missing_dark_constants} are not available for correction.")
if corr_bools.get('relgain') and "RelativeGainCCD" not in constants.keys():
warning("RelativeGainEPix100 is not found in the calibration database.")
corr_bools['relgain'] = False
step_timer.done_step("Constants retrieval")
```
%% Cell type:code id: tags:
``` python
fig = xana.heatmapPlot(constants["Offset"][:,:,0], x_label='Columns', y_label='Rows', lut_label='Offset (ADU)',
fig = xana.heatmapPlot(constants["OffsetCCD"][:,:,0], x_label='Columns', y_label='Rows', lut_label='Offset (ADU)',
aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x), vmax=16000,
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Dark Offset Map')
fig = xana.heatmapPlot(constants["Noise"][:,:,0], x_label='Columns', y_label='Rows',
fig = xana.heatmapPlot(constants["NoiseCCD"][:,:,0], x_label='Columns', y_label='Rows',
lut_label='Corrected Noise (ADU)',
aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Dark Noise Map')
fig = xana.heatmapPlot(np.log2(constants["BadPixelsDark"][:,:,0]), x_label='Columns', y_label='Rows',
fig = xana.heatmapPlot(np.log2(constants["BadPixelsDarkCCD"][:,:,0]), x_label='Columns', y_label='Rows',
lut_label='Bad Pixel Value (ADU)',
aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Dark Bad Pixels Map')
if corr_bools.get('relgain'):
fig = xana.heatmapPlot(constants["RelativeGain"], figsize=(8, 8), x_label='Columns', y_label='Rows',
fig = xana.heatmapPlot(constants["RelativeGainCCD"], figsize=(8, 8), x_label='Columns', y_label='Rows',
lut_label='Relative Gain',
aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=0.8, vmax=1.2,
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
panel_top_low_lim = 0.5, panel_top_high_lim = 1.5, panel_side_low_lim = 0.5,
panel_side_high_lim = 1.5,
title = f'Relative Gain Map for pnCCD (Gain = 1/{int(gain)})')
step_timer.done_step("Constants retrieval")
```
%% Cell type:code id: tags:
``` python
#************************ Calculators ************************#
if corr_bools.get('common_mode'):
# Common Mode Correction Calculator:
cmCorrection = xcal.CommonModeCorrection([pixels_x, pixels_y],
commonModeBlockSize,
commonModeAxis,
parallel=False, dType=np.float32, stride=1,
noiseMap=constants["Noise"].astype(np.float32), minFrac=0.25)
noiseMap=constants["NoiseCCD"].astype(np.float32), minFrac=0.25)
if corr_bools.get('pattern_class'):
# Pattern Classifier Calculator:
# Left Hemisphere:
patternClassifierLH = xcal.PatternClassifier([pixels_x, pixels_y//2],
constants["Noise"][:, :pixels_y//2],
constants["NoiseCCD"][:, :pixels_y//2],
split_evt_primary_threshold,
split_evt_secondary_threshold,
split_evt_mip_threshold,
tagFirstSingles=3, # track along y-axis, left to right (see
nCells=1, # split_event.py file in pydetlib/lib/src/
allowElongated=False, # XFELDetAna/algorithms)
blockSize=[pixels_x, pixels_y//2],
parallel=False)
# Right Hemisphere:
patternClassifierRH = xcal.PatternClassifier([pixels_x, pixels_y//2],
constants["Noise"][:, pixels_y//2:],
constants["NoiseCCD"][:, pixels_y//2:],
split_evt_primary_threshold,
split_evt_secondary_threshold,
split_evt_mip_threshold,
tagFirstSingles=4, # track along y-axis, right to left
nCells=1,
allowElongated=False,
blockSize=[pixels_x, pixels_y//2],
parallel=False)
patternClassifierLH._imagesPerChunk = 1
patternClassifierRH._imagesPerChunk = 1
patternClassifierLH._noisemap = constants["Noise"][:, :pixels_x//2]
patternClassifierRH._noisemap = constants["Noise"][:, pixels_x//2:]
patternClassifierLH._noisemap = constants["NoiseCCD"][:, :pixels_x//2]
patternClassifierRH._noisemap = constants["NoiseCCD"][:, pixels_x//2:]
# Setting bad pixels:
patternClassifierLH.setBadPixelMask(constants["BadPixelsDark"][:, :pixels_x//2] != 0)
patternClassifierRH.setBadPixelMask(constants["BadPixelsDark"][:, pixels_x//2:] != 0)
patternClassifierLH.setBadPixelMask(constants["BadPixelsDarkCCD"][:, :pixels_x//2] != 0)
patternClassifierRH.setBadPixelMask(constants["BadPixelsDarkCCD"][:, pixels_x//2:] != 0)
```
%% Cell type:code id: tags:
``` python
#***************** Histogram Calculators ******************#
# Will contain uncorrected data:
histCalRaw = xcal.HistogramCalculator(sensorSize,
bins=bins,
range=bin_range,
nCells=1,
parallel=False,
blockSize=blockSize)
# Will contain offset corrected data:
histCalOffsetCor = xcal.HistogramCalculator(sensorSize,
bins=bins,
range=bin_range,
nCells=1,
parallel=False,
blockSize=blockSize)
if corr_bools.get('common_mode'):
# Will contain common mode corrected data:
histCalCommonModeCor = xcal.HistogramCalculator(sensorSize,
bins=bins,
range=bin_range,
nCells=1,
parallel=False,
blockSize=blockSize)
if corr_bools.get('pattern_class'):
# Will contain split events pattern data:
histCalPcorr = xcal.HistogramCalculator(sensorSize,
bins=bins,
range=bin_range,
nCells=1,
parallel=False,
blockSize=blockSize)
# Will contain singles events data:
histCalPcorrS = xcal.HistogramCalculator(sensorSize,
bins=bins,
range=bin_range,
nCells=1,
parallel=False,
blockSize=blockSize)
if corr_bools.get('relgain'):
# Will contain gain corrected data:
histCalGainCor = xcal.HistogramCalculator(sensorSize,
bins=bins,
range=bin_range,
nCells=1,
parallel=False,
blockSize=blockSize)
```
%% Cell type:markdown id: tags:
## Applying corrections to the raw data
%% Cell type:code id: tags:
``` python
def offset_correction(wid, index, d):
"""offset correction.
Equating bad pixels' values to np.nan,
so that the pattern classifier ignores them:
"""
d = d.copy()
# TODO: To clear this up. Is it on purpose to save corrected data with nans?
d[bpix != 0] = np.nan
d -= offset # offset correction
# TODO: to clear this up. why save the badpixels map in the corrected data?
bpix_data[index, ...] = bpix
data[index, ...] = d
def common_mode(wid, index, d):
"""common-mode correction.
Discarding events caused by saturated pixels:
"""
d = np.squeeze(cmCorrection.correct(d, cellTable=np.zeros(pixels_y, np.int32)))
# we equate these values to np.nan so that the pattern classifier ignores them:
d[d >= saturated_threshold] = np.nan
data[index, ...] = d
def gain_correction(wid, index, d):
"""relative gain correction."""
d /= relativegain
data[index, ...] = d
def pattern_classification_correction(wid, index, d):
"""pattern classification correction.
data set to save split event corrected images
The calculation of the cluster map:]
Dividing the data into left and right hemispheres:
"""
# pattern classification on corrected data
dataLH, patternsLH = patternClassifierLH.classify(d[:, :pixels_x//2])
dataRH, patternsRH = patternClassifierRH.classify(d[:, pixels_x//2:])
d[:, :pixels_x//2] = np.squeeze(dataLH)
d[:, pixels_x//2:] = np.squeeze(dataRH)
patterns = np.zeros(d.shape, patternsLH.dtype)
patterns[:, :pixels_x//2] = np.squeeze(patternsLH)
patterns[:, pixels_x//2:] = np.squeeze(patternsRH)
d[d < split_evt_primary_threshold*noise] = 0
data[index, ...] = d
ptrn_data[index, ...] = patterns
d[patterns != 100] = np.nan # Discard doubles, triples, quadruple, clusters, first singles
filtered_data[index, ...] = d
```
%% Cell type:code id: tags:
``` python
# 10 is a number chosen after testing 1 ... 71 parallel threads for a node with 72 cpus.
parallel_num_threads = 10
context = psh.context.ThreadContext(num_workers=parallel_num_threads)
data_path = "INSTRUMENT/"+instrument_src+"/data/"
offset = np.squeeze(constants["Offset"])
noise = np.squeeze(constants["Noise"])
bpix = np.squeeze(constants["BadPixelsDark"])
relativegain = constants.get("RelativeGain")
offset = np.squeeze(constants["OffsetCCD"])
noise = np.squeeze(constants["NoiseCCD"])
bpix = np.squeeze(constants["BadPixelsDarkCCD"])
relativegain = constants.get("RelativeGainCCD")
```
%% Cell type:code id: tags:
``` python
def write_datasets(corr_arrays, ofile):
"""
Creating datasets first then adding data.
To have metadata together available at the start of the file,
so it's quick to see what the file contains
"""
comp_fields = ["gain", "patterns", "pixels_classified"]
img_grp = ofile[data_path]
for field, data in corr_arrays.items():
kw = dict(chunks=(chunk_size_idim, pixels_x, pixels_y))
if field in comp_fields:
kw["compression"] = "gzip"
img_grp.create_dataset(
field, shape=data.shape, dtype=data.dtype, **kw)
for field, data in corr_arrays.items():
img_grp[field][:] = data
```
%% Cell type:code id: tags:
``` python
# Data corrections and event classifications happen here.
# Also, the corrected data are written to datasets:
for seq_n, seq_f in enumerate(seq_files):
f_dc = H5File(seq_f)
out_file = f"{out_folder}/{seq_f.name}".replace("RAW", "CORR")
step_timer.start()
dshape = f_dc[instrument_src, "data.image"].shape
n_trains = dshape[0]
# If you want to analyze only a certain number of frames
# instead of all available good frames.
if limit_images > 0:
n_trains = min(n_trains, limit_images)
data_shape = (n_trains, dshape[1], dshape[2])
print(f"Correcting file: {seq_f} of shape {data_shape}.")
data_dc = f_dc.select(instrument_src, "data.image", require_all=True).select_trains(np.s_[:n_trains]) # noqa
raw_data = data_dc[instrument_src, "data.image"].ndarray().astype(np.float32)
if seq_n == 0:
raw_plt = raw_data.copy() # plot first sequence only
step_timer.start()
# Allocating shared arrays for data arrays for each correction stage.
data = context.alloc(shape=data_shape, dtype=np.float32)
bpix_data = context.alloc(shape=data_shape, dtype=np.uint32)
histCalRaw.fill(raw_data) # filling histogram with raw uncorrected data
# Applying offset correction
context.map(offset_correction, raw_data)
histCalOffsetCor.fill(data) # filling histogram with offset corrected data
if seq_n == 0:
off_data = data.copy() # plot first sequence only
corr_arrays = {
"pixels": data.copy(),
"mask": bpix_data,
}
step_timer.done_step(f'offset correction.')
if corr_bools.get('common_mode'):
step_timer.start()
# Applying common mode correction
context.map(common_mode, data)
if seq_n == 0:
cm_data = data.copy() # plot first sequence only
corr_arrays["pixels_cm"] = data.copy()
histCalCommonModeCor.fill(data) # filling histogram with common mode corrected data
step_timer.done_step(f'common-mode correction.')
if corr_bools.get('relgain'):
step_timer.start()
# Applying gain correction
context.map(gain_correction, data)
if seq_n == 0:
rg_data = data.copy() # plot first sequence only
# TODO: Why storing a repeated constant for each image in corrected files.
corr_arrays["gain"] = np.repeat(relativegain[np.newaxis, ...], n_trains, axis=0).astype(np.float32) # noqa
histCalGainCor.fill(data) # filling histogram with gain corrected data
step_timer.done_step(f'gain correction.')
if corr_bools.get('pattern_class'):
step_timer.start()
ptrn_data = context.alloc(shape=data_shape, dtype=np.int32)
filtered_data = context.alloc(shape=data_shape, dtype=np.int32)
# Applying pattern classification correction
# Even thougth data is indeed of dtype np.float32,
# not specifiying this again screw with the data quality.
context.map(pattern_classification_correction, data.astype(np.float32))
if seq_n == 0:
cls_data = data.copy() # plot first sequence only
# split event corrected images plotted for first sequence only
# (also these events are only singles events):
corr_arrays["pixels_classified"] = data.copy()
corr_arrays["patterns"] = ptrn_data
histCalPcorr.fill(data) # filling histogram with split events corrected data
# filling histogram with corr data after discarding doubles, triples, quadruple, clusters, and first singles
histCalPcorrS.fill(filtered_data)
step_timer.done_step(f'pattern classification correction.')
step_timer.start()
# Storing corrected data sources.
with h5py.File(out_file, 'w') as ofile:
# Copy RAW non-calibrated sources.
with h5py.File(seq_f, 'r') as sfile:
h5_copy_except.h5_copy_except_paths(
sfile, ofile, [],
)
# TODO: to clear this up: why save corrected data in data/pixels rather than data/image.
write_datasets(corr_arrays, ofile)
step_timer.done_step(f'Storing data.')
```
%% Cell type:code id: tags:
``` python
print("In addition to offset correction, the following corrections were performed:")
for k, v in corr_bools.items():
if v:
print(" -", k.upper())
print(f"Total processing time {step_timer.timespan():.01f} s")
step_timer.print_summary()
```
%% Cell type:code id: tags:
``` python
print("In addition to offset correction, the following corrections were performed:")
for k, v in corr_bools.items():
if v:
print(" -", k.upper())
print(f"Total processing time {step_timer.timespan():.01f} s")
step_timer.print_summary()
```
%% Cell type:code id: tags:
``` python
# Histograming the resulting spectra:
# First _ refers to the bin edges and second _ refers to statistics and we ignore both.
# if you use histCalRaw.get(cumulatative = True) and so on, the cumulatative = True turns the counts array such as
# RawHistVals and so on into a 1D array instead of keeping the original shape:
RawHistVals, _, RawHistMids, _ = histCalRaw.get()
off_cor_HistVals, _, off_cor_HistMids, _ = histCalOffsetCor.get()
if corr_bools.get('common_mode'):
cm_cor_HistVals, _, cm_HistMids, _ = histCalCommonModeCor.get()
if corr_bools.get('relgain'):
gain_cor_HistVals, _, gain_cor_HistMids, _ = histCalGainCor.get()
if corr_bools.get('pattern_class'):
split_HistVals, _, split_HistMids, _ = histCalPcorr.get() # split events corrected
singles_HistVals, _, singles_HistMids, _ = histCalPcorrS.get() # last s in variable names: singles events
```
%% Cell type:code id: tags:
``` python
# Saving intermediate data to disk:
step_timer.start()
np.savez(os.path.join(out_folder, 'Raw_Events.npz'), RawHistMids, RawHistVals)
np.savez(os.path.join(out_folder, 'Offset_Corrected_Events.npz'), off_cor_HistMids, off_cor_HistVals)
if corr_bools.get('common_mode'):
np.savez(os.path.join(out_folder, 'Common_Mode_Corrected_Events.npz'), cm_HistMids, cm_cor_HistVals)
if corr_bools.get('relgain'):
np.savez(os.path.join(out_folder, 'Gain_Corrected_Events.npz'), gain_cor_HistMids, gain_cor_HistVals)
if corr_bools.get('pattern_class'):
np.savez(os.path.join(out_folder, 'Split_Events_Corrected_Events.npz'), split_HistMids, split_HistVals)
np.savez(os.path.join(out_folder, 'Singles_Events.npz'), singles_HistMids, singles_HistVals)
step_timer.done_step(f'Saving intermediate data to disk.')
print("Various spectra are saved to disk in the form of histograms. Please check {}".format(out_folder))
```
%% Cell type:code id: tags:
``` python
display(Markdown('### Raw vs. Corrected Spectra'))
step_timer.start()
figure = [{'x': RawHistMids,
'y': RawHistVals,
'y_err': np.sqrt(RawHistVals[:]),
'drawstyle': 'steps-post',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'Uncorrected'
},
{'x': off_cor_HistMids,
'y': off_cor_HistVals,
'y_err': np.sqrt(off_cor_HistVals[:]),
'drawstyle': 'steps-post',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'Offset Corrected'
}]
if corr_bools.get('common_mode'):
figure.append({'x': cm_HistMids,
'y': cm_cor_HistVals,
'y_err': np.sqrt(cm_cor_HistVals[:]),
'drawstyle': 'steps-post',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'Common Mode Corrected'})
if corr_bools.get('relgain'):
xrange = bin_range
figure.append({'x': gain_cor_HistMids,
'y': gain_cor_HistVals,
'y_err': np.sqrt(gain_cor_HistVals[:]),
'drawstyle': 'steps-post',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'Gain Corrected'})
if corr_bools.get('pattern_class'):
figure.extend([{'x': split_HistMids,
'y': split_HistVals,
'y_err': np.sqrt(split_HistVals[:]),
'drawstyle': 'steps-post',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'Split Events Corrected'
},
{'x': singles_HistMids,
'y': singles_HistVals,
'y_err': np.sqrt(singles_HistVals[:]),
'drawstyle': 'steps-post',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'Singles Events'
}])
fig = xana.simplePlot(figure, aspect=1, x_label='ADU', y_label='Number of Occurrences', figsize='2col',
y_log=True, x_range=bin_range, title = '1 ADU per bin is used.',
legend='top-right-frame-1col')
step_timer.done_step('Plotting')
```
%% Cell type:code id: tags:
``` python
# This function plots pattern statistics:
def classification_plot(patternStats, hemisphere):
print("****************** {} HEMISPHERE ******************\n"
.format(hemisphere))
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(4, 4, 1)
sfields = ["singles", "first singles", "clusters"]
mfields = ["doubles", "triples", "quads"]
relativeOccurances = []
labels = []
for i, f in enumerate(sfields):
relativeOccurances.append(patternStats[f])
labels.append(f)
for i, f in enumerate(mfields):
for k in range(len(patternStats[f])):
relativeOccurances.append(patternStats[f][k])
labels.append("{}({})".format(f, k))
relativeOccurances = np.array(relativeOccurances, np.float)
relativeOccurances /= np.sum(relativeOccurances)
pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)
ax.set_title("Pattern Occurrence")
# Set aspect ratio to be equal so that pie is drawn as a circle.
a = ax.axis('equal')
smaps = ["singlemap", "firstsinglemap", "clustermap"]
for i, m in enumerate(smaps):
ax = fig.add_subplot(4, 4, 2+i)
pmap = ax.imshow(patternStats[m], interpolation="nearest", vmax=2*np.nanmedian(patternStats[m]))
ax.set_title(m)
cb = fig.colorbar(pmap)
mmaps = ["doublemap", "triplemap", "quadmap"]
k = 0
for i, m in enumerate(mmaps):
for j in range(4):
ax = fig.add_subplot(4, 4, 2+len(smaps)+k)
pmap = ax.imshow(patternStats[m][j], interpolation="nearest", vmax=2*np.median(patternStats[m][j]))
ax.set_title("{}({})".format(m,j))
cb = fig.colorbar(pmap)
k+=1
```
%% Cell type:code id: tags:
``` python
# The next two cells plot the classification results for left and right hemispheres, respectively:
display(Markdown('### Classification Results - Plots'))
if corr_bools.get('pattern_class'):
patternStatsLH = patternClassifierLH.getPatternStats()
classification_plot(patternStatsLH, 'Left')
```
%% Cell type:code id: tags:
``` python
if corr_bools.get('pattern_class'):
patternStatsRH = patternClassifierRH.getPatternStats()
classification_plot(patternStatsRH, 'Right')
```
%% Cell type:code id: tags:
``` python
display(Markdown('### Classification Results - Tabulated Statistics'))
if corr_bools.get('pattern_class'):
step_timer.start()
t0 = PrettyTable()
t0.title = "Total Number of Counts after All Corrections"
t0.field_names = ["Hemisphere", "Singles", "First-Singles", "Clusters"]
t0.add_row(["LH", patternStatsLH['singles'], patternStatsLH['first singles'], patternStatsLH['clusters']])
t0.add_row(["RH", patternStatsRH['singles'], patternStatsRH['first singles'], patternStatsRH['clusters']])
print(t0)
print("Abbreviations: D (Doubles), T (Triples), Q (Quadruples), L (Left), R (Right), and H (Hemisphere).")
t1 = PrettyTable()
t1.field_names = ["Index", "D-LH", "D-RH", "T-LH", "T-RH", "Q-LH", "Q-RH"]
t1.add_row([0, patternStatsLH['doubles'][0], patternStatsRH['doubles'][0], patternStatsLH['triples'][0],
patternStatsRH['triples'][0], patternStatsLH['quads'][0], patternStatsRH['quads'][0]])
t1.add_row([1, patternStatsLH['doubles'][1], patternStatsRH['doubles'][1], patternStatsLH['triples'][1],
patternStatsRH['triples'][1], patternStatsLH['quads'][1], patternStatsRH['quads'][1]])
t1.add_row([2, patternStatsLH['doubles'][2], patternStatsRH['doubles'][2], patternStatsLH['triples'][2],
patternStatsRH['triples'][2], patternStatsLH['quads'][2], patternStatsRH['quads'][2]])
t1.add_row([3, patternStatsLH['doubles'][3], patternStatsRH['doubles'][3], patternStatsLH['triples'][3],
patternStatsRH['triples'][3], patternStatsLH['quads'][3], patternStatsRH['quads'][3]])
print(t1)
step_timer.done_step('Classification Results - Tabulated Statistics')
```
%% Cell type:code id: tags:
``` python
if corr_bools.get('pattern_class'):
doublesLH = patternStatsLH['doubles'][0] + patternStatsLH['doubles'][1] + patternStatsLH['doubles'][2] + \
patternStatsLH['doubles'][3]
triplesLH = patternStatsLH['triples'][0] + patternStatsLH['triples'][1] + patternStatsLH['triples'][2] + \
patternStatsLH['triples'][3]
quadsLH = patternStatsLH['quads'][0] + patternStatsLH['quads'][1] + patternStatsLH['quads'][2] + \
patternStatsLH['quads'][3]
allsinglesLH = patternStatsLH['singles'] + patternStatsLH['first singles']
eventsLH = allsinglesLH + doublesLH + triplesLH + quadsLH
doublesRH = patternStatsRH['doubles'][0] + patternStatsRH['doubles'][1] + patternStatsRH['doubles'][2] + \
patternStatsRH['doubles'][3]
triplesRH = patternStatsRH['triples'][0] + patternStatsRH['triples'][1] + patternStatsRH['triples'][2] + \
patternStatsRH['triples'][3]
quadsRH = patternStatsRH['quads'][0] + patternStatsRH['quads'][1] + patternStatsRH['quads'][2] + \
patternStatsRH['quads'][3]
allsinglesRH = patternStatsRH['singles'] + patternStatsRH['first singles']
eventsRH = allsinglesRH + doublesRH + triplesRH + quadsRH
if eventsLH > 0.:
reloccurLH = np.array([allsinglesLH/eventsLH, doublesLH/eventsLH, triplesLH/eventsLH, quadsLH/eventsLH])
else:
reloccurLH = np.array([0]*4)
if eventsRH > 0.:
reloccurRH = np.array([allsinglesRH/eventsRH, doublesRH/eventsRH, triplesRH/eventsRH, quadsRH/eventsRH])
else:
reloccurRH = np.array([0]*4)
```
%% Cell type:code id: tags:
``` python
display(Markdown('### Classification Results - Pie Charts'))
if corr_bools.get('pattern_class'):
step_timer.start()
fig = plt.figure(figsize=(12, 7))
ax = fig.add_subplot(1, 2, 1)
labels = ['Singles', 'Doubles', 'Triples', 'Quads']
pie = ax.pie(reloccurLH, labels=labels, autopct='%1.1f%%', shadow=True)
ax.set_title("Pattern Occurrence in LH")
# Set aspect ratio to be equal so that pie is drawn as a circle.
a = ax.axis('equal')
ax = fig.add_subplot(1, 2, 2)
pie = ax.pie(reloccurRH, labels=labels, autopct='%1.1f%%', shadow=True)
ax.set_title("Pattern Occurrence in RH")
# Set aspect ratio to be equal so that pie is drawn as a circle.
a = ax.axis('equal')
step_timer.done_step('Classification Results - Pie Charts')
```
%% Cell type:markdown id: tags:
### Various Images Averaged Over All Frames of Only the First Sequence ###
%% Cell type:code id: tags:
``` python
step_timer.start()
uncor_mean_im = np.nanmean(raw_data, axis=0)
offset_mean_im = np.nanmean(off_data, axis=0)
if corr_bools.get('common_mode'):
cm_mean_im = np.nanmean(cm_data, axis=0)
if corr_bools.get('relgain'):
gain_mean_im = np.nanmean(rg_data, axis=0)
if corr_bools.get('pattern_class'):
mean_im_cc = np.nanmean(cls_data, axis=0)
fig = xana.heatmapPlot(uncor_mean_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Uncorrected Image Averaged over Frames in the First Sequence')
fig = xana.heatmapPlot(offset_mean_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Offset Corrected Image Averaged over Frames in the First Sequence')
if corr_bools.get('common_mode'):
fig = xana.heatmapPlot(cm_mean_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Common Mode Corrected Image Averaged over Frames in the First Sequence')
if corr_bools.get('relgain'):
fig = xana.heatmapPlot(gain_mean_im, x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)', aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Gain Corrected Image Averaged over Frames in the First Sequence')
if corr_bools.get('pattern_class'):
fig = xana.heatmapPlot(mean_im_cc, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=0, vmax= 18000,
title = 'Image of Single Events Averaged over Frames in the First Sequence')
step_timer.done_step("Plotting")
```
%% Cell type:markdown id: tags:
### Images of the First Frame of the First Sequence ###
%% Cell type:code id: tags:
``` python
step_timer.start()
fig = xana.heatmapPlot(raw_data[0, :, :], x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Uncorrected Image (First Frame of the First Sequence)')
fig = xana.heatmapPlot(off_data[0, :, :], x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Offset Corrected Image (First Frame of the First Sequence)')
if corr_bools.get('common_mode'):
fig = xana.heatmapPlot(cm_data[0, :, :], x_label='Columns', y_label='Rows', lut_label='Signal (ADU)',
aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Common Mode Corrected Image (First Frame of the First Sequence)')
if corr_bools.get('relgain'):
fig = xana.heatmapPlot(rg_data[0, :, :], x_label='Columns', y_label='Rows', lut_label='Signal (ADU)',
aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Gain Corrected Image (First Frame of the First Sequence)')
if corr_bools.get('pattern_class'):
fig = xana.heatmapPlot(cls_data[0, :, :], x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
x_range=(0, pixels_y), y_range=(0, pixels_x),
panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
title = 'Image of Single Events (First Frame of the First Sequence)')
step_timer.done_step("Plotting")
```
%% Cell type:code id: tags:
``` python
# Resetting the histogram calculators:
histCalRaw.reset()
histCalOffsetCor.reset()
if corr_bools.get('common_mode'):
histCalCommonModeCor.reset()
if corr_bools.get('relgain'):
histCalGainCor.reset()
if corr_bools.get('pattern_class'):
histCalPcorr.reset()
histCalPcorrS.reset()
```
%% Cell type:markdown id: tags:
Next, the corrected event patterns are read from the patterns/dataset created previously and are separated into 4 different categories (singles, doubles, triples and quadruples) using the pattern indices. However, this is done only for one sequence, corresponding to the seq_num variable, as an example.
Note that the number of bins and the bin range for the following histograms may be different from those presented above (depending on gain) to make the counts more noticible and the peaks more defined.
If you are interested in plotting the events from all sequences or the spectra of half of the sensor, execute the spectra_pnCCD_NBC.ipynb notebook.
%% Cell type:code id: tags:
``` python
if corr_bools.get('pattern_class'):
singles = []
doubles = []
triples = []
quads = []
with H5File(f"{out_folder}/{seq_files[0].name.replace('RAW', 'CORR')}") as dc: # noqa
data = dc[instrument_src, "data.pixels_classified"].ndarray()
patterns = dc[instrument_src, "data.patterns"].ndarray()
# events' patterns indices are as follows: 100 (singles), 101 (first singles), 200 - 203 (doubles),
# 300 - 303 (triples), and 400 - 403 (quadruples). Note that for the last three types of patterns,
# there are left, right, up, and down indices.
# Separating the events:
# Singles and First Singles:
for s in range(100, 102):
single = data.copy()
single[patterns != s] = np.nan
singles.append(single)
for d in range(200, 204):
double = data.copy()
double[patterns != d] = np.nan
doubles.append(double)
for t in range(300, 304):
triple = data.copy()
triple[patterns != t] = np.nan
triples.append(triple)
for q in range(400, 404):
quad = data.copy()
quad[patterns != q] = np.nan
quads.append(quad)
```
%% Cell type:code id: tags:
``` python
if corr_bools.get('pattern_class'):
step_timer.start()
hA = 0
h = 0
for single in singles:
hs, e = np.histogram(single.flatten(), bins=event_bins, range=b_range) # h: histogram counts, e: bin edges
h += hs
hA += hs # hA: counts all events (see below)
# bin edges array has one extra element => need to plot from 0 to the one before the last element to have the
# same size as h-array => in what follows, we use e[:-1] (-1 means one before the last element)
display(Markdown('### Histograms of Corrected Events for One Sequence Only'))
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111)
ax.step(e[:-1], h, color='blue', label='Events Involving Single Pixels Only')
ax.semilogy() # y-axis is log, x-axis is linear
ax.set_xlabel("Energy (ADU) [{} bins per {} ADU]".format(event_bins, b_range[1]-b_range[0]))
ax.set_ylabel("Corrected Events for One Sequence (counts)")
ax.set_xlim(b_range)
h = 0
for double in doubles:
hd, e = np.histogram(double.flatten(), bins=event_bins, range=b_range)
h += hd
hA += hd
ax.step(e[:-1], h, color='red', label='Events Splitting on Double Pixels')
h = 0
for triple in triples:
ht, e = np.histogram(triple.flatten(), bins=event_bins, range=b_range)
h += ht
hA += ht
ax.step(e[:-1], h, color='green', label='Events Splitting on Triple Pixels')
h = 0
for quad in quads:
hq, e = np.histogram(quad.flatten(), bins=event_bins, range=b_range)
h += hq
hA += hq
ax.step(e[:-1], h, color='purple', label='Events Splitting on Quadruple Pixels')
ax.step(e[:-1], hA, color='grey', label='All Valid Events')
l = ax.legend()
step_timer.done_step("Plotting")
```
%% Cell type:code id: tags:
``` python
print(f"Total processing time {step_timer.timespan():.01f} s")
step_timer.print_summary()
```

Legal Notice | Contacts