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

Check if Relgain constant exists

parent 03244729
No related branches found
No related tags found
1 merge request!203Feat&fix: (FastCCD CORRECTION) add bools to control correction steps
Hide whitespace changes
Inline Side-by-side
notebooks/FastCCD/CorrectionNotebook_NewDAQ_FastCCD_NBC.ipynb
+ 18
− 14
View file @ 8c8aae1f
%% Cell type:markdown id: tags:
# FastCCD Data Correction ##
Authors: I. Klačková, S. Hauf, Version 1.0
The following notebook provides correction of images acquired with the FastCCD.
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/DETLAB/201931/p900104/raw/" # input folder, required
out_folder = '/gpfs/exfel/data/scratch/ahmedk/test/fastccd' # output folder, required
in_folder = "/gpfs/exfel/exp/SCS/201930/p900074/raw" # input folder, required
out_folder = '/gpfs/exfel/data/scratch/karnem/test/fastccd' # output folder, required
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # path template in hdf5 file
path_inset = 'DA01'
run = 114 # run number
path_inset = 'DA05'
run = 453 # run number
h5path = '/INSTRUMENT/SCS_CDIDET_FCCD2M/DAQ/FCCD:daqOutput/data/image' # path in HDF5 file
h5path_t = '/CONTROL/SCS_CDIDET_FCCD2M/CTRL/LSLAN/inputA/crdg/value' # temperature path in HDF5 file
h5path_cntrl = '/RUN/SCS_CDIDET_FCCD2M/DET/FCCD' # path to control data
cluster_profile = "noDB" #ipcluster profile to use
cpuCores = 16 #Specifies the number of running cpu cores
operation_mode = "FF" # FS stands for frame-store and FF for full-frame opeartion
split_evt_primary_threshold = 7. # primary threshold for split event classification in terms of n sigma noise
split_evt_secondary_threshold = 4. # secondary threshold for split event classification in terms of n sigma noise
split_evt_mip_threshold = 1000. # MIP threshold for event classification
cal_db_interface = "tcp://max-exfl016:8015#8025" # calibration DB interface to use
cal_db_timeout = 300000000 # timeout on caldb requests
sequences = [-1] # sequences to correct, set to -1 for all, range allowed
chunk_size_idim = 1 # H5 chunking size of output data
overwrite = True # overwrite existing files
sequences_per_node = 1 # sequences to correct per node
limit_images = 0 # limit images per file
use_dir_creation_date = True # use dir creation data for calDB queries
time_offset_days = 0 # offset in days for calibration parameters
photon_energy_gain_map = 5.9 # energy in keV
fix_temperature = 0. # fix temperature to this value, set to 0 to use slow control value
flipped_between = ["2019-02-01", "2019-04-02"] # detector was flipped during this timespan
temp_limits = 5 # limits within which temperature is considered the same
commonModeAxis = 1 # Axis along which common mode will be calculated (0: along rows, 1: along columns)
# Correction Booleans
only_offset = False # Only apply offset correction
cti_corr = False # Apply CTI correction
relgain_corr = True # Apply relative gain correction
common_mode_corr = False # Apply commonMode correction
correct_offset_drift = False # correct for offset drifts
do_pattern_classification = True # classify split events
def balance_sequences(in_folder, run, sequences, sequences_per_node, path_inset):
import glob
import re
import numpy as np
if sequences[0] == -1:
sequence_files = glob.glob("{}/r{:04d}/*{}-S*.h5".format(in_folder, run, path_inset))
seq_nums = set()
for sf in sequence_files:
seqnum = re.findall(r".*-S([0-9]*).h5", sf)[0]
seq_nums.add(int(seqnum))
seq_nums -= set(sequences)
nsplits = len(seq_nums)//sequences_per_node+1
while nsplits > 8:
sequences_per_node += 1
nsplits = len(seq_nums)//sequences_per_node+1
print("Changed to {} sequences per node to have a maximum of 8 concurrent jobs".format(sequences_per_node))
return [l.tolist() for l in np.array_split(list(seq_nums), nsplits)]
else:
return sequences
```
%% Cell type:code id: tags:
``` python
# Fill dictionaries comprising bools and arguments for correction and data analysis
# Here the herarichy and dependability for correction booleans are defined
corr_bools = {}
# offset is at the bottom of AGIPD correction pyramid.
corr_bools["only_offset"] = only_offset
# Dont apply any corrections if only_offset is requested
if not only_offset:
# Apply relative gain correction, only if requested
if relgain_corr:
corr_bools["relgain"] = relgain_corr
# Apply CTI correction, only if requested
if cti_corr:
corr_bools["cti"] = cti_corr
corr_bools["common_mode"] = common_mode_corr
corr_bools["offset_drift"] = correct_offset_drift
corr_bools["pattern_class"] = do_pattern_classification
# Here the herarichy and dependability for data analysis booleans and arguments are defined
data_analysis_parms = {}
```
%% Cell type:code id: tags:
``` python
import XFELDetAna.xfelprofiler as xprof
profiler = xprof.Profiler()
profiler.disable()
from XFELDetAna.util import env
env.iprofile = cluster_profile
import warnings
warnings.filterwarnings('ignore')
from XFELDetAna import xfelpycaltools as xcal
from XFELDetAna import xfelpyanatools as xana
from XFELDetAna.plotting.util import prettyPlotting
prettyPlotting=True
from XFELDetAna.xfelreaders import ChunkReader
from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5
import numpy as np
import h5py
import matplotlib.pyplot as plt
from iminuit import Minuit
import time
import copy
import os
from prettytable import PrettyTable
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from iCalibrationDB.detectors import DetectorTypes
from cal_tools.tools import (get_dir_creation_date,
get_random_db_interface,
get_constant_from_db_and_time)
from datetime import timedelta
%matplotlib inline
if sequences[0] == -1:
sequences = None
# select a random port for the data base
cal_db_interface = get_random_db_interface(cal_db_interface)
```
%% Cell type:code id: tags:
``` python
if operation_mode == "FS":
x = 960 # rows of the FastCCD to analyze in FS mode
y = 960 # columns of the FastCCD to analyze in FS mode
print('\nYou are analyzing data in FS mode.')
else:
x = 1934 # rows of the FastCCD to analyze in FF mode
y = 960 # columns of the FastCCD to analyze in FF mode
print('\nYou are analyzing data in FF mode.')
ped_dir = "{}/r{:04d}".format(in_folder, run)
out_folder = "{}/r{:04d}".format(out_folder, run)
fp_name = path_template.format(run, path_inset)
fp_path = '{}/{}'.format(ped_dir, fp_name)
print("Reading data from: {}\n".format(fp_path))
print("Run is: {}".format(run))
print("HDF5 path: {}".format(h5path))
print("Data is output to: {}".format(out_folder))
import datetime
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run) + timedelta(days=time_offset_days)
if creation_time:
print("Using {} as creation time".format(creation_time.isoformat()))
```
%% Cell type:code id: tags:
``` python
sensorSize = [x, y]
chunkSize = 100 #Number of images to read per chunk
blockSize = [sensorSize[0]//2, sensorSize[1]] #Sensor area will be analysed according to blocksize
xcal.defaultBlockSize = blockSize
memoryCells = 1 #FastCCD has 1 memory cell
#Specifies total number of images to proceed
commonModeBlockSize = blockSize
# commonModeAxisR = 'row'#Axis along which common mode will be calculated
run_parallel = True
profile = False
filename = fp_path.format(sequences[0] if sequences else 0)
with h5py.File(filename, 'r') as f:
bias_voltage = int(f['{}/biasclock/bias/value'.format(h5path_cntrl)][0])
det_gain = int(f['{}/exposure/gain/value'.format(h5path_cntrl)][0])
integration_time = int(f['{}/exposure/exposure_time/value'.format(h5path_cntrl)][0])
print("Bias voltage is {} V".format(bias_voltage))
print("Detector gain is set to x{}".format(det_gain))
print("Detector integration time is set to {}".format(integration_time))
temperature = np.mean(f[h5path_t])
temperature_k = temperature + 273.15
if fix_temperature != 0.:
temperature_k = fix_temperature
print("Using fixed temperature")
print("Mean temperature was {:0.2f} °C / {:0.2f} K at beginning of run".format(temperature, temperature_k))
if not os.path.exists(out_folder):
os.makedirs(out_folder)
elif not overwrite:
# Stop Notebook not only this cell
raise SystemExit("Output path exists! Exiting")
```
%% Cell type:code id: tags:
``` python
dirlist = sorted(os.listdir(ped_dir))
file_list = []
total_sequences = 0
fsequences = []
for entry in dirlist:
#only h5 file
abs_entry = "{}/{}".format(ped_dir, entry)
if os.path.isfile(abs_entry) and os.path.splitext(abs_entry)[1] == ".h5":
if sequences is None:
for seq in range(len(dirlist)):
if path_template.format(run, path_inset).format(seq) in abs_entry:
file_list.append(abs_entry)
total_sequences += 1
fsequences.append(seq)
else:
for seq in sequences:
if path_template.format(run, path_inset).format(seq) in abs_entry:
file_list.append(os.path.abspath(abs_entry))
total_sequences += 1
fsequences.append(seq)
sequences = fsequences
```
%% Cell type:code id: tags:
``` python
import copy
from IPython.display import HTML, display, Markdown, Latex
import tabulate
print("Processing a total of {} sequence files".format(total_sequences))
table = []
for k, f in enumerate(file_list):
table.append((k, f))
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
offsetMap = None
badPixelMap = None
noiseMap = None
# The following constants are saved in data base as a constant for each gain
# Hence, a loop over all 3 gains is performed
# This is to be used in messages in reports.
g_name = {8: "High gain", 2: "Medium gain", 1: "Low gain"}
for i, g in enumerate([8, 2, 1]):
print("Retrieving constants for {}".format(g_name[g]))
# set the operating condition
condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
integration_time=integration_time,
gain_setting=g,
temperature=temperature_k,
pixels_x=x,
pixels_y=y)
for parm in condition.parameters:
if parm.name == "Sensor Temperature":
parm.lower_deviation = temp_limits
parm.upper_deviation = temp_limits
offset = Constants.CCD(DetectorTypes.fastCCD).Offset()
noise = Constants.CCD(DetectorTypes.fastCCD).Noise()
bpix = Constants.CCD(DetectorTypes.fastCCD).BadPixelsDark()
## retrieve_offset
offset, offset_time = get_constant_from_db_and_time(device=Detectors.fastCCD1,
constant=offset,
condition=condition,
empty_constant=None,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout, print_once=False)
if offsetMap is None:
offsetMap = np.zeros(list(offset.shape)+[3], np.float32)
if offset is not None:
offsetMap[...,i] = offset
else:
print("NO OFFSET FOUND IN DB!")
offset_temperature = None
for parm in condition.parameters:
if parm.name == "Sensor Temperature":
offset_temperature = parm.value
print("Dark Offset was taken at temperature of {:0.2f}K at {}"
.format(offset_temperature, offset_time))
## retrieve_noise
noise, noise_time = get_constant_from_db_and_time(device=Detectors.fastCCD1,
constant=noise,
condition=condition,
empty_constant=None,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout, print_once=False)
if noiseMap is None:
noiseMap = np.zeros(list(noise.shape)+[3], np.float32)
if noise is not None:
noiseMap[...,i] = noise
else:
print("NO NOISE FOUND IN DB!")
print("Noise at {} was taken at {}"
.format(g_name[g], noise_time))
## retrieve_bad pixels
bpix, bpix_time = get_constant_from_db_and_time(device=Detectors.fastCCD1,
constant=bpix,
condition=condition,
empty_constant=None,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout, print_once=False)
if badPixelMap is None:
badPixelMap = np.zeros(list(bpix.shape)+[3], np.uint32)
print("NO BadPixel FOUND IN DB!")
if bpix is not None:
badPixelMap[...,i] = noise
else:
print("NO BADPIXEL FOUND IN DB!")
print("BadPixes at {} was taken at {}"
.format(g_name[g], bpix_time))
```
%% Cell type:markdown id: tags:
Loading cti and relative gain values
%% Cell type:code id: tags:
``` python
## relative gain
# relative gain
if corr_bools.get('relgain'):
relgain = Constants.CCD(DetectorTypes.fastCCD).RelativeGain()
# set the operating condition
condition = Conditions.Illuminated.CCD(bias_voltage=bias_voltage,
integration_time=integration_time,
gain_setting=det_gain,
temperature=temperature_k,
pixels_x=x, pixels_y=y,
photon_energy=photon_energy_gain_map)
relgain, relgain_time = get_constant_from_db_and_time(device=Detectors.fastCCD1,
constant=relgain,
condition=condition,
empty_constant=None,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout, print_once=False)
constant=relgain,
condition=condition,
empty_constant=None,
cal_db_interface=cal_db_interface,
creation_time=creation_time,
timeout=cal_db_timeout, print_once=False)
# TODO: CHECK IF THIS FLIPPING IS CORRECT
relGain = relgain[::-1,...]
print('Relative Gain was taken with creation-date:', relgain_time)
if relgain is None:
corr_bools["relgain"] = False
print('Relative Gain was not found. Proceed without Relative Gain correction.')
else:
relGain = relgain[::-1, ...]
print('Relative Gain was taken with creation-date:', relgain_time)
```
%% Cell type:code id: tags:
``` python
if corr_bools.get('cti'):
pass
## FASTCCD CTI CURRENTLY IS NOT SUPPORTED BY DET!
# The code is left for tracing past algorithm for later
# use during the development of the new CTI algorithm.
# TODO: Proper CTI Retrieval from CTI and the correction
# in following cells
# relGainCA = copy.copy(relGain)
# relGainC = relGainCA[:relGainCA.shape[0]//2,...]
# ctiA = np.ones(relGainCA.shape[:2])
# cti = np.ones(relGainC.shape[:2])
# i = 0
# idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
# mn1 = np.nanmean(relGainC[i, ~idx, 0])
# for i in range(1, relGainC.shape[0]):
# idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
# mn2 = np.nanmean(relGainC[i, ~idx, 0])
# cti[i,:] = mn2/mn1
# ctiA[:relGainCA.shape[0]//2,...] = cti
# relGainC = relGainCA[relGainCA.shape[0]//2:,...]
# cti = np.ones(relGainC.shape[:2])
# i = -1
# idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
# mn1 = np.nanmean(relGainC[i, ~idx, 0])
# for i in range(relGainC.shape[0]-1, 1, -1):
# idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
# mn2 = np.nanmean(relGainC[i, ~idx, 0])
# cti[i,:] = mn2/mn1
# ctiA[relGainCA.shape[0]//2:,...] = cti
# relGainCA = copy.copy(relGain)
# relGainC = relGainCA[:relGainCA.shape[0]//2,...]
# for i in range(relGainC.shape[1]):
# idx = (relGainC[:,i, 0] < 0.95) | (relGainC[:,i,0] > 1.05)
# relGainC[idx,i,0] = np.nanmean(relGainC[~idx,i,0])
# relGainC[idx,i,1] = np.nanmean(relGainC[~idx,i,1])
# relGainC[idx,i,2] = np.nanmean(relGainC[~idx,i,2])
# relGainCA[:relGainCA.shape[0]//2,...] = relGainC
# relGainC = relGainCA[relGainCA.shape[0]//2:,...]
# for i in range(relGainC.shape[1]):
# idx = (relGainC[:,i, 0] < 0.95) | (relGainC[:,i,0] > 1.05)
# relGainC[idx,i,0] = np.nanmean(relGainC[~idx,i,0])
# relGainC[idx,i,1] = np.nanmean(relGainC[~idx,i,1])
# relGainC[idx,i,2] = np.nanmean(relGainC[~idx,i,2])
# relGainCA[relGainCA.shape[0]//2:,...] = relGainC
# relGainC = relGainCA*ctiA[...,None]
# relGain = relGainC
```
%% Cell type:code id: tags:
``` python
import dateutil.parser
flipped_between = [dateutil.parser.parse(d) for d in flipped_between]
flip_rgain = creation_time.replace(tzinfo=None) >= flipped_between[0] and creation_time.replace(tzinfo=None) <= flipped_between[1]
flip_rgain &= (relgain_time.replace(tzinfo=None) >= flipped_between[0]
and relgain_time.replace(tzinfo=None) <= flipped_between[1])
print("Accounting for flipped detector: {}".format(flip_rgain))
```
%% Cell type:code id: tags:
``` python
#************************Calculators************************#
if corr_bools.get('common_mode'):
cmCorrection = xcal.CommonModeCorrection([x, y],
commonModeBlockSize,
commonModeAxis,
nCells = memoryCells,
stride=10,
runParallel=True,
stats=True, minFrac=0)
patternClassifierLH = xcal.PatternClassifier([x//2, y],
noiseMap[:x//2, :],
split_evt_primary_threshold,
split_evt_secondary_threshold,
split_evt_mip_threshold,
tagFirstSingles = 0,
nCells=memoryCells,
cores=cpuCores,
allowElongated = False,
blockSize=[x//2, y],
runParallel=True)
patternClassifierUH = xcal.PatternClassifier([x//2, y],
noiseMap[x//2:, :],
split_evt_primary_threshold,
split_evt_secondary_threshold,
split_evt_mip_threshold,
tagFirstSingles = 0,
nCells=memoryCells,
cores=cpuCores,
allowElongated = False,
blockSize=[x//2, y],
runParallel=True)
```
%% Cell type:code id: tags:
``` python
#*****************Histogram Calculators******************#
histCalOffsetCor = xcal.HistogramCalculator([x, y],
bins=1050,
range=[-50, 1000],
nCells=memoryCells,
cores=cpuCores,
blockSize=blockSize)
histCalPcorr = xcal.HistogramCalculator([x, y],
bins=1050,
range=[-50, 1000],
nCells=memoryCells,
cores=cpuCores,
blockSize=blockSize)
histCalPcorrS = xcal.HistogramCalculator([x, y],
bins=1050,
range=[-50, 1000],
nCells=memoryCells,
cores=cpuCores,
blockSize=blockSize)
histCalCommonModeCor = xcal.HistogramCalculator([x, y],
bins=1050,
range=[-50, 1000],
nCells=memoryCells,
cores=cpuCores,
blockSize=blockSize)
```
%% Cell type:markdown id: tags:
Applying corrections
%% Cell type:code id: tags:
``` python
patternClassifierLH._imagesPerChunk = 500
patternClassifierUH._imagesPerChunk = 500
patternClassifierLH.debug()
patternClassifierUH.debug()
```
%% Cell type:code id: tags:
``` python
histCalOffsetCor.debug()
histCalCommonModeCor.debug()
histCalPcorr.debug()
```
%% Cell type:code id: tags:
``` python
def copy_and_sanitize_non_cal_data(infile, outfile, h5base):
"""
:param infile: Input file
:param outfile: Otput file
:param h5base:
"""
if h5base.startswith("/"):
h5base = h5base[1:]
dont_copy = ['pixels']
dont_copy = [h5base+"/{}".format(do)
for do in dont_copy]
def visitor(k, item):
"""
:param k:
:param 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
mean_im = None
single_im = None
mean_im_cc = None
single_im_cc = None
drift_lh = []
drift_uh = []
offsetMap = np.squeeze(offsetMap)
noiseMap = np.squeeze(noiseMap)
badPixelMap = np.squeeze(badPixelMap)
if corr_bools.get('relgain'):
#TODO: This should be removed after properly injecting gain const for all 3 gains
if len(relGain.shape) == 3 and relGain.shape[2] == 1:
# This is a temporary solution for the relGain of one gain in db
relGain = np.repeat(relGain, 3, axis=2)
relGain = np.squeeze(relGain)
for k, f in enumerate(file_list):
with h5py.File(f, 'r', driver='core') as infile:
out_fileb = "{}/{}".format(out_folder, f.split("/")[-1])
out_file = out_fileb.replace("RAW", "CORR")
#out_filed = out_fileb.replace("RAW", "CORR-SC")
data = None
noise = None
try:
with h5py.File(out_file, "w") as ofile:
copy_and_sanitize_non_cal_data(infile, ofile, h5path)
data = infile[h5path+"/pixels"][()]
nzidx = np.count_nonzero(data, axis=(1,2))
data = data[nzidx != 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)
ddsetm = ofile.create_dataset(h5path+"/mask",
oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.uint32, compression="gzip")
ddsetg = ofile.create_dataset(h5path+"/gain",
oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.uint8, compression="gzip")
# Getting the 14th and 15th bit from the data,
# which contains gain before removing them from the data.
gain = np.right_shift(data, 14)
# gains are stored as 00 for High gain,
# 01 for Medium gain and 11 for low gain.
# Hence, the subtraction from the
# gain's int values to have 0, 1 and 2
gain[gain != 0] -= 1
fstride = 1
if not flip_rgain: # rgain was taken during flipped orientation
fstride = -1
data = np.bitwise_and(data, 0b0011111111111111).astype(np.float32)
# creating maps for correction usage
omap = np.repeat(offsetMap[...,None,:], data.shape[2], axis=2)
nmap = np.repeat(noiseMap[...,None,:], data.shape[2], axis=2)
bmap = np.repeat(badPixelMap[...,None,:], data.shape[2], axis=2)
# selecting element value related to the gain of the pixel.
offset = np.choose(gain, (omap[...,0], omap[...,1], omap[...,2]))
noise = np.choose(gain, (nmap[...,0], nmap[...,1], nmap[...,2]))
bpix = np.choose(gain, (bmap[...,0], bmap[...,1], bmap[...,2]))
# same for relative gain if correction is available
if corr_bools.get('relgain'):
rmap = np.repeat(relGain[:,::fstride,None,:], data.shape[2], axis=2)
rg = np.choose(gain, (rmap[...,0], rmap[...,1], rmap[...,2]))
# Apply offset correction
data -= offset
if corr_bools.get('relgain'):
# Apply relative gain correction
# TODO: check relgain correction in pydetlib
# and use it here if the same.
data *= rg
if corr_bools.get("offset_drift"):
# Put in consideration the temperature
# change of the FastCCD's hole.
lhd = np.mean(data[x//2-10:x//2,y//2-5:y//2+5,:], axis=(0,1))
data[:x//2, :, :] -= lhd
drift_lh.append(lhd)
uhd = np.mean(data[x//2:x//2+10,y//2-5:y//2+5,:], axis=(0,1))
data[x//2:, :, :] -= uhd
drift_uh.append(uhd)
histCalOffsetCor.fill(data)
ddset[...] = np.moveaxis(data, 2, 0)
ddsetm[...] = np.moveaxis(bpix, 2, 0)
ddsetg[...] = np.moveaxis(gain, 2, 0).astype(np.uint8)
if mean_im is None:
mean_im = np.nanmean(data, axis=2)
single_im = data[...,0]
if corr_bools.get("pattern_class"):
ddsetcm = ofile.create_dataset(h5path+"/pixels_cm",
oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.float32)
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")
# The calculation of the cluster map
patternClassifierLH._noisemap = noise[:x//2, :, :]
patternClassifierUH._noisemap = noise[x//2:, :, :]
# common mode correction
if corr_bools.get("common_mode"):
cellTable = np.zeros(data.shape[2], np.int32) # Common mode correction
data = cmCorrection.correct(data.astype(np.float32),
cellTable=cellTable,
noiseMap=noise)
# correct for the row common mode
ddsetcm[...] = np.moveaxis(data, 2, 0)
histCalCommonModeCor.fill(data)
dataLH = data[:x//2, :, :]
dataUH = data[x//2:, :, :]
dataLH, patternsLH = patternClassifierLH.classify(dataLH)
dataUH, patternsUH = patternClassifierUH.classify(dataUH)
data[:x//2, :, :] = dataLH
data[x//2:, :, :] = dataUH
patterns = np.zeros(data.shape, patternsLH.dtype)
patterns[:x//2, :, :] = patternsLH
patterns[x//2:, :, :] = patternsUH
data[data < split_evt_primary_threshold*noise] = 0
ddsetc[...] = np.moveaxis(data, 2, 0)
ddsetp[...] = np.moveaxis(patterns, 2, 0)
histCalPcorr.fill(data)
data[patterns != 100] = np.nan
histCalPcorrS.fill(data)
if mean_im_cc is None:
mean_im_cc = np.nanmean(data, axis=2)
single_im_cc = data[...,0]
except Exception as e:
print("Couldn't calibrate data in {}: {}".format(f, e))
```
%% Cell type:code id: tags:
``` python
if corr_bools.get("offset_drift"):
lhds = np.concatenate(drift_lh)
uhds = np.concatenate(drift_uh)
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(111)
ax.plot(lhds, label="Lower hem.")
ax.plot(uhds, label="Upper hem.")
ax.set_xlabel("Frame #")
ax.set_xlabel("Offset drift (ADU)")
```
%% Cell type:code id: tags:
``` python
if corr_bools.get("pattern_class"):
print("******************LOWER HEMISPHERE******************\n")
patternStatsLH = patternClassifierLH.getPatternStats()
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(patternStatsLH[f])
labels.append(f)
for i, f in enumerate(mfields):
for k in range(len(patternStatsLH[f])):
relativeOccurances.append(patternStatsLH[f][k])
labels.append(f+"("+str(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(patternStatsLH[m], interpolation="nearest", vmax=2*np.nanmedian(patternStatsLH[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(patternStatsLH[m][j], interpolation="nearest", vmax=2*np.median(patternStatsLH[m][j]))
ax.set_title(m+"("+str(j)+")")
cb = fig.colorbar(pmap)
k+=1
```
%% Cell type:code id: tags:
``` python
if corr_bools.get("pattern_class"):
print("******************UPPER HEMISPHERE******************\n")
patternStatsUH = patternClassifierUH.getPatternStats()
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(patternStatsUH[f])
labels.append(f)
for i, f in enumerate(mfields):
for k in range(len(patternStatsUH[f])):
relativeOccurances.append(patternStatsUH[f][k])
labels.append(f+"("+str(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(patternStatsUH[m], interpolation="nearest", vmax=2*np.nanmedian(patternStatsUH[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(patternStatsUH[m][j], interpolation="nearest", vmax=np.median(patternStatsUH[m][j]))
ax.set_title(m+"("+str(j)+")")
cb = fig.colorbar(pmap)
k+=1
```
%% Cell type:code id: tags:
``` python
if corr_bools.get("pattern_class"):
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(["UH", patternStatsUH['singles'], patternStatsUH['first singles'], patternStatsUH['clusters']])
print(t0)
t1 = PrettyTable()
t1.field_names = ["Index","D-LH", "D-UH", "T-LH", "T-UH", "Q-LH", "Q-UH"]
t1.add_row([0, patternStatsLH['doubles'][0], patternStatsUH['doubles'][0], patternStatsLH['triples'][0], patternStatsUH['triples'][0], patternStatsLH['quads'][0], patternStatsUH['quads'][0]])
t1.add_row([1, patternStatsLH['doubles'][1], patternStatsUH['doubles'][1], patternStatsLH['triples'][1], patternStatsUH['triples'][1], patternStatsLH['quads'][1], patternStatsUH['quads'][1]])
t1.add_row([2, patternStatsLH['doubles'][2], patternStatsUH['doubles'][2], patternStatsLH['triples'][2], patternStatsUH['triples'][2], patternStatsLH['quads'][2], patternStatsUH['quads'][2]])
t1.add_row([3, patternStatsLH['doubles'][3], patternStatsUH['doubles'][3], patternStatsLH['triples'][3], patternStatsUH['triples'][3], patternStatsLH['quads'][3], patternStatsUH['quads'][3]])
print(t1)
```
%% 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
doublesUH = patternStatsUH['doubles'][0] + patternStatsUH['doubles'][1] + patternStatsUH['doubles'][2] + patternStatsUH['doubles'][3]
triplesUH = patternStatsUH['triples'][0] + patternStatsUH['triples'][1] + patternStatsUH['triples'][2] + patternStatsUH['triples'][3]
quadsUH = patternStatsUH['quads'][0] + patternStatsUH['quads'][1] + patternStatsUH['quads'][2] + patternStatsUH['quads'][3]
allsinglesUH = patternStatsUH['singles'] + patternStatsUH['first singles']
eventsUH = allsinglesUH + doublesUH + triplesUH + quadsUH
reloccurLH = np.array([allsinglesLH/eventsLH, doublesLH/eventsLH, triplesLH/eventsLH, quadsLH/eventsLH])
reloccurUH = np.array([allsinglesUH/eventsUH, doublesUH/eventsUH, triplesUH/eventsUH, quadsUH/eventsUH])
```
%% Cell type:code id: tags:
``` python
if corr_bools.get("pattern_class"):
fig = plt.figure(figsize=(10,5))
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 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(reloccurUH, labels=labels, autopct='%1.1f%%', shadow=True)
ax.set_title("Pattern occurrence UH")
# Set aspect ratio to be equal so that pie is drawn as a circle.
a = ax.axis('equal')
```
%% 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 corr_bools.get("pattern_class") and corr_bools.get("common_mode"):
hcm,ecm,ccm,scm = histCalCommonModeCor.get()
d.append({'x': ccm,
'y': hcm,
'y_err': np.sqrt(hcm[:]),
'drawstyle': 'steps-mid',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'CommonMode corr.'
})
fig = xana.simplePlot(d, aspect=1, x_label='Energy(ADU)',
y_label='Number of occurrences', figsize='2col',
y_log=True, x_range=(-50,500),
legend='top-center-frame-2col')
```
%% Cell type:code id: tags:
``` python
if corr_bools.get("pattern_class"):
h1,e1L,c1L,s1L = histCalPcorr.get()
h1s,e1Ls,c1Ls,s1Ls = histCalPcorrS.get()
d = [
{'x': c1L,
'y': h1,
'y_err': np.sqrt(h1[:]),
'drawstyle': 'steps-mid',
'label': 'Split event corrected'},
{'x': c1Ls,
'y': h1s,
'y_err': np.sqrt(h1s[:]),
'drawstyle': 'steps-mid',
'label': 'Single pixel hits'}
]
fig = xana.simplePlot(d, aspect=1, x_label='Energy(ADU)',
y_label='Number of occurrences', figsize='2col',
y_log=True, x_range=(0,200),x_log=False,
legend='top-center-frame-2col')
```
%% Cell type:markdown id: tags:
## Mean Image of first Sequence ##
%% Cell type:code id: tags:
``` python
fig = xana.heatmapPlot(mean_im,
x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)',
x_range=(0,y),
y_range=(0,x), vmin=-50, vmax=500)
if corr_bools.get("pattern_class"):
fig = xana.heatmapPlot(mean_im_cc,
x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)',
x_range=(0,y),
y_range=(0,x), vmin=-50, vmax=500)
```
%% Cell type:markdown id: tags:
## Single Shot of first Sequnce ##
%% Cell type:code id: tags:
``` python
fig = xana.heatmapPlot(single_im,
x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)',
x_range=(0,y),
y_range=(0,x), vmin=-50, vmax=500)
if corr_bools.get("pattern_class"):
fig = xana.heatmapPlot(single_im_cc,
x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)',
x_range=(0,y),
y_range=(0,x), vmin=-50, vmax=500)
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Legal Notice | Contacts