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

Merge branch 'fix/cm_correction' into 'master' 

Fix: PNCCD DARK, CM correction

See merge request detectors/pycalibration!239
parents 409fb398 a4a65246
No related branches found
No related tags found
1 merge request!239Fix: PNCCD DARK, CM correction
Hide whitespace changes
Inline Side-by-side
notebooks/pnCCD/Correct_pnCCD_NBC.ipynb
+ 43
− 65
View file @ 97a5afe3
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# pnCCD Data Correction ## # pnCCD Data Correction ##
Authors: DET Group, Version 1.0 Authors: DET Group, Version 1.0
The following notebook provides correction of images acquired with the pnCCD. The following notebook provides correction of images acquired with the pnCCD.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
in_folder = "/gpfs/exfel/exp/SQS/201921/p002430/raw/" # input folder, required in_folder = "/gpfs/exfel/exp/SQS/201921/p002430/raw/" # input folder, required
out_folder = '/gpfs/exfel/data/scratch/xcal/test/' # output folder, required out_folder = '/gpfs/exfel/data/scratch/karnem/test/' # output folder, required
path_template = 'RAW-R{:04d}-PNCCD01-S{{:05d}}.h5' # the template to use to access data path_template = 'RAW-R{:04d}-PNCCD01-S{{:05d}}.h5' # the template to use to access data
path_template_seqs = "{}/r{:04d}/*PNCCD01-S*.h5" path_template_seqs = "{}/r{:04d}/*PNCCD01-S*.h5"
run = 746 # which run to read data from, required run = 746 # which run to read data from, required
number_dark_frames = 0 # number of images to be used, if set to 0 all available images are used number_dark_frames = 0 # number of images to be used, if set to 0 all available images are used
cluster_profile = "noDB" # ipcluster profile to use cluster_profile = "noDB" # ipcluster profile to use
sigma_noise = 10. # Pixel exceeding 'sigmaNoise' * noise value in that pixel will be masked sigma_noise = 10. # Pixel exceeding 'sigmaNoise' * noise value in that pixel will be masked
h5path = '/INSTRUMENT/SQS_NQS_PNCCD1MP/CAL/PNCCD_FMT-0:output/data/' # path in the HDF5 file the data is at h5path = '/INSTRUMENT/SQS_NQS_PNCCD1MP/CAL/PNCCD_FMT-0:output/data/' # path in the HDF5 file the data is at
multi_iteration = False # use multiple iterations multi_iteration = False # use multiple iterations
use_dir_creation_date = True use_dir_creation_date = True
fix_temperature = 233. fix_temperature = 233.
gain = 0 gain = 0
bias_voltage = 300 bias_voltage = 300
integration_time = 70 integration_time = 70
split_evt_primary_threshold = 5. # primary threshold for split event classification in terms of n sigma noise split_evt_primary_threshold = 5. # 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 split_evt_secondary_threshold = 3. # secondary threshold for split event classification in terms of n sigma noise
split_evt_mip_threshold = 1000. # MIP threshold for event classification split_evt_mip_threshold = 1000. # MIP threshold for event classification
cal_db_interface = "tcp://max-exfl016:8015" # calibration DB interface to use cal_db_interface = "tcp://max-exfl016:8015" # calibration DB interface to use
cal_db_timeout = 300000000 # timeout on caldb requests cal_db_timeout = 300000000 # timeout on caldb requests
sequences = [-1] # sequences to correct, set to -1 for all, range allowed sequences = [-1] # sequences to correct, set to -1 for all, range allowed
chunk_size_idim = 1 # H5 chunking size of output data chunk_size_idim = 1 # H5 chunking size of output data
overwrite = True overwrite = True
do_pattern_classification = True # classify split events do_pattern_classification = True # classify split events
sequences_per_node = 1 sequences_per_node = 1
limit_images = 0 limit_images = 0
time_offset_days = 0 time_offset_days = 0
photon_energy_gain_map = 2. # energy in keV photon_energy_gain_map = 2. # energy in keV
cpuCores = 8 cpuCores = 8
def balance_sequences(in_folder, run, sequences, sequences_per_node, path_template_seqs): def balance_sequences(in_folder, run, sequences, sequences_per_node, path_template_seqs):
import glob import glob
import re import re
import numpy as np import numpy as np
if sequences[0] == -1: if sequences[0] == -1:
sequence_files = glob.glob(path_template_seqs.format(in_folder, run)) sequence_files = glob.glob(path_template_seqs.format(in_folder, run))
seq_nums = set() seq_nums = set()
for sf in sequence_files: for sf in sequence_files:
seqnum = re.findall(r".*-S([0-9]*).h5", sf)[0] seqnum = re.findall(r".*-S([0-9]*).h5", sf)[0]
seq_nums.add(int(seqnum)) seq_nums.add(int(seqnum))
seq_nums -= set(sequences) seq_nums -= set(sequences)
nsplits = len(seq_nums)//sequences_per_node+1 nsplits = len(seq_nums)//sequences_per_node+1
while nsplits > 8: while nsplits > 8:
sequences_per_node += 1 sequences_per_node += 1
nsplits = len(seq_nums)//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)) 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)] return [l.tolist() for l in np.array_split(list(seq_nums), nsplits)]
else: else:
return sequences return sequences
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import XFELDetAna.xfelprofiler as xprof import XFELDetAna.xfelprofiler as xprof
profiler = xprof.Profiler() profiler = xprof.Profiler()
profiler.disable() profiler.disable()
from XFELDetAna.util import env from XFELDetAna.util import env
env.iprofile = cluster_profile env.iprofile = cluster_profile
import warnings import warnings
warnings.filterwarnings('ignore') warnings.filterwarnings('ignore')
from XFELDetAna import xfelpycaltools as xcal from XFELDetAna import xfelpycaltools as xcal
from XFELDetAna import xfelpyanatools as xana from XFELDetAna import xfelpyanatools as xana
from XFELDetAna.plotting.util import prettyPlotting from XFELDetAna.plotting.util import prettyPlotting
prettyPlotting=True prettyPlotting=True
from XFELDetAna.xfelreaders import ChunkReader from XFELDetAna.xfelreaders import ChunkReader
from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5 from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5
import numpy as np import numpy as np
import h5py import h5py
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from iminuit import Minuit from iminuit import Minuit
import time import time
import copy import copy
import os import os
from prettytable import PrettyTable from prettytable import PrettyTable
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from iCalibrationDB.detectors import DetectorTypes from iCalibrationDB.detectors import DetectorTypes
from cal_tools.tools import get_dir_creation_date from cal_tools.tools import get_dir_creation_date
from datetime import timedelta from datetime import timedelta
%matplotlib inline %matplotlib inline
if sequences[0] == -1: if sequences[0] == -1:
sequences = None sequences = None
if "#" in cal_db_interface: if "#" in cal_db_interface:
prot, serv, ran = cal_db_interface.split(":") prot, serv, ran = cal_db_interface.split(":")
r1, r2 = ran.split("#") r1, r2 = ran.split("#")
cal_db_interface = ":".join( cal_db_interface = ":".join(
[prot, serv, str(np.random.randint(int(r1), int(r2)))]) [prot, serv, str(np.random.randint(int(r1), int(r2)))])
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
x = 1024 # rows of the FastCCD to analyze in FS mode x = 1024 # rows of the FastCCD to analyze in FS mode
y = 1024 # columns of the FastCCD to analyze in FS mode y = 1024 # columns of the FastCCD to analyze in FS mode
ped_dir = "{}/r{:04d}".format(in_folder, run) ped_dir = "{}/r{:04d}".format(in_folder, run)
fp_name = path_template.format(run) fp_name = path_template.format(run)
import datetime import datetime
creation_time = None creation_time = None
if use_dir_creation_date: if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run) creation_time = get_dir_creation_date(in_folder, run)
fp_path = '{}/{}'.format(ped_dir, fp_name) fp_path = '{}/{}'.format(ped_dir, fp_name)
print("Reading data from: {}\n".format(fp_path)) print("Reading data from: {}\n".format(fp_path))
print("Run is: {}".format(run)) print("Run is: {}".format(run))
print("HDF5 path: {}".format(h5path)) print("HDF5 path: {}".format(h5path))
if creation_time: if creation_time:
print("Using {} as creation time".format(creation_time.isoformat())) print("Using {} as creation time".format(creation_time.isoformat()))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
sensorSize = [x, y] sensorSize = [x, y]
chunkSize = 100 #Number of images to read per chunk chunkSize = 100 #Number of images to read per chunk
blockSize = [sensorSize[0]//2, sensorSize[1]//4] #Sensor area will be analysed according to blocksize blockSize = [sensorSize[0]//2, sensorSize[1]//4] #Sensor area will be analysed according to blocksize
xcal.defaultBlockSize = blockSize xcal.defaultBlockSize = blockSize
memoryCells = 1 #FastCCD has 1 memory cell memoryCells = 1 #FastCCD has 1 memory cell
#Specifies total number of images to proceed #Specifies total number of images to proceed
commonModeBlockSize = [512, 512] commonModeBlockSize = [512, 512]
commonModeAxisR = 'row'#Axis along which common mode will be calculated commonModeAxisR = 'row'#Axis along which common mode will be calculated
run_parallel = True run_parallel = True
profile = False profile = False
if not os.path.exists(out_folder): if not os.path.exists(out_folder):
os.makedirs(out_folder) os.makedirs(out_folder)
elif not overwrite: elif not overwrite:
raise AttributeError("Output path exists! Exiting") raise AttributeError("Output path exists! Exiting")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
dirlist = sorted(os.listdir(ped_dir)) dirlist = sorted(os.listdir(ped_dir))
file_list = [] file_list = []
total_sequences = 0 total_sequences = 0
fsequences = [] fsequences = []
for entry in dirlist: for entry in dirlist:
#only h5 file #only h5 file
abs_entry = "{}/{}".format(ped_dir, entry) abs_entry = "{}/{}".format(ped_dir, entry)
if os.path.isfile(abs_entry) and os.path.splitext(abs_entry)[1] == ".h5": if os.path.isfile(abs_entry) and os.path.splitext(abs_entry)[1] == ".h5":
if sequences is None: if sequences is None:
for seq in range(len(dirlist)): for seq in range(len(dirlist)):
if path_template.format(run).format(seq) in abs_entry: if path_template.format(run).format(seq) in abs_entry:
file_list.append(abs_entry) file_list.append(abs_entry)
total_sequences += 1 total_sequences += 1
fsequences.append(seq) fsequences.append(seq)
else: else:
for seq in sequences: for seq in sequences:
if path_template.format(run).format(seq) in abs_entry: if path_template.format(run).format(seq) in abs_entry:
file_list.append(os.path.abspath(abs_entry)) file_list.append(os.path.abspath(abs_entry))
total_sequences += 1 total_sequences += 1
fsequences.append(seq) fsequences.append(seq)
sequences = fsequences sequences = fsequences
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import copy import copy
from IPython.display import HTML, display, Markdown, Latex from IPython.display import HTML, display, Markdown, Latex
import tabulate import tabulate
print("Processing a total of {} sequence files".format(total_sequences)) print("Processing a total of {} sequence files".format(total_sequences))
table = [] table = []
for k, f in enumerate(file_list): for k, f in enumerate(file_list):
table.append((k, f)) table.append((k, f))
if len(table): if len(table):
md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=["#", "file"]))) md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=["#", "file"])))
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
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
offsetMap = None offsetMap = None
offset_temperature = None offset_temperature = None
badPixelMap = None badPixelMap = None
noiseMap = None noiseMap = None
## offset ## offset
metadata = ConstantMetaData() metadata = ConstantMetaData()
offset = Constants.CCD(DetectorTypes.pnCCD).Offset() offset = Constants.CCD(DetectorTypes.pnCCD).Offset()
metadata.calibration_constant = offset metadata.calibration_constant = offset
# set the operating condition # set the operating condition
condition = Conditions.Dark.CCD(bias_voltage=bias_voltage, condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
integration_time=integration_time, integration_time=integration_time,
gain_setting=gain, gain_setting=gain,
temperature=fix_temperature, temperature=fix_temperature,
pixels_x=1024, pixels_x=1024,
pixels_y=1024) pixels_y=1024)
device = Detectors.PnCCD1 device = Detectors.PnCCD1
metadata.detector_condition = condition metadata.detector_condition = condition
# specify the a version for this constant # specify the a version for this constant
if creation_time is None: if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device) metadata.calibration_constant_version = Versions.Now(device=device)
metadata.retrieve(cal_db_interface) metadata.retrieve(cal_db_interface)
else: else:
metadata.calibration_constant_version = Versions.Timespan(device=device, metadata.calibration_constant_version = Versions.Timespan(device=device,
start=creation_time) start=creation_time)
metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000) metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)
if offsetMap is None: if offsetMap is None:
offsetMap = np.zeros(list(offset.data.shape)+[3], np.float32) offsetMap = np.zeros(list(offset.data.shape)+[3], np.float32)
offsetMap = offset.data offsetMap = offset.data
offset_temperature = None offset_temperature = None
for parm in condition.parameters: for parm in condition.parameters:
if parm.name == "Sensor Temperature": if parm.name == "Sensor Temperature":
offset_temperature = parm.value offset_temperature = parm.value
print("Temperature dark images for offset calculation " + print("Temperature dark images for offset calculation " +
"were taken at: {:0.2f} K @ {}".format(offset_temperature, "were taken at: {:0.2f} K @ {}".format(offset_temperature,
metadata.calibration_constant_version.begin_at)) metadata.calibration_constant_version.begin_at))
## noise ## noise
metadata = ConstantMetaData() metadata = ConstantMetaData()
noise = Constants.CCD(DetectorTypes.pnCCD).Noise() noise = Constants.CCD(DetectorTypes.pnCCD).Noise()
metadata.calibration_constant = noise metadata.calibration_constant = noise
# set the operating condition # set the operating condition
condition = Conditions.Dark.CCD(bias_voltage=bias_voltage, condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
integration_time=integration_time, integration_time=integration_time,
gain_setting=gain, gain_setting=gain,
temperature=fix_temperature, temperature=fix_temperature,
pixels_x=1024, pixels_x=1024,
pixels_y=1024) pixels_y=1024)
device = Detectors.PnCCD1 device = Detectors.PnCCD1
metadata.detector_condition = condition metadata.detector_condition = condition
# specify the a version for this constant # specify the a version for this constant
if creation_time is None: if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device) metadata.calibration_constant_version = Versions.Now(device=device)
metadata.retrieve(cal_db_interface) metadata.retrieve(cal_db_interface)
else: else:
metadata.calibration_constant_version = Versions.Timespan(device=device, metadata.calibration_constant_version = Versions.Timespan(device=device,
start=creation_time) start=creation_time)
metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000) metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)
noiseMap = noise.data noiseMap = noise.data
## bad pixels ## bad pixels
metadata = ConstantMetaData() metadata = ConstantMetaData()
bpix = Constants.CCD(DetectorTypes.pnCCD).BadPixelsDark() bpix = Constants.CCD(DetectorTypes.pnCCD).BadPixelsDark()
metadata.calibration_constant = bpix metadata.calibration_constant = bpix
# set the operating condition # set the operating condition
condition = Conditions.Dark.CCD(bias_voltage=bias_voltage, condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
integration_time=integration_time, integration_time=integration_time,
gain_setting=gain, gain_setting=gain,
temperature=fix_temperature, temperature=fix_temperature,
pixels_x=1024, pixels_x=1024,
pixels_y=1024) pixels_y=1024)
device = Detectors.PnCCD1 device = Detectors.PnCCD1
metadata.detector_condition = condition metadata.detector_condition = condition
# specify the a version for this constant # specify the a version for this constant
metadata.calibration_constant_version = Versions.Now(device=device) metadata.calibration_constant_version = Versions.Now(device=device)
metadata.retrieve(cal_db_interface, timeout=cal_db_timeout) metadata.retrieve(cal_db_interface, timeout=cal_db_timeout)
badPixelMap = bpix.data badPixelMap = bpix.data
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Loading cti and relative gain values Loading cti and relative gain values
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
## relative gain ## relative gain
metadata = ConstantMetaData() metadata = ConstantMetaData()
relgain = Constants.CCD(DetectorTypes.pnCCD).RelativeGain() relgain = Constants.CCD(DetectorTypes.pnCCD).RelativeGain()
metadata.calibration_constant = relgain metadata.calibration_constant = relgain
# set the operating condition # set the operating condition
condition = Conditions.Illuminated.CCD(bias_voltage=bias_voltage, condition = Conditions.Illuminated.CCD(bias_voltage=bias_voltage,
integration_time=integration_time, integration_time=integration_time,
gain_setting=gain, gain_setting=gain,
temperature=fix_temperature, temperature=fix_temperature,
pixels_x=1024, pixels_x=1024,
pixels_y=1024, photon_energy=photon_energy_gain_map) pixels_y=1024, photon_energy=photon_energy_gain_map)
device = Detectors.PnCCD1 device = Detectors.PnCCD1
metadata.detector_condition = condition metadata.detector_condition = condition
# specify the a version for this constant # specify the a version for this constant
metadata.calibration_constant_version = Versions.Now(device=device) metadata.calibration_constant_version = Versions.Now(device=device)
#metadata.retrieve(cal_db_interface) #metadata.retrieve(cal_db_interface)
relGain = np.ones_like(offsetMap) #relgain.data relGain = np.ones_like(offsetMap) #relgain.data
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#************************Calculators************************# #************************Calculators************************#
cmCorrection = xcal.CommonModeCorrection([x, y], cmCorrection = xcal.CommonModeCorrection([x, y],
commonModeBlockSize, commonModeBlockSize,
commonModeAxisR, commonModeAxisR,
nCells = memoryCells, parallel=False, dType=np.float32, stride=1,
noiseMap = noiseMap, noiseMap=noiseMap.astype(np.float32), minFrac=0)
runParallel=True, # cmCorrection.debug()
stats=True)
#cmCorrection.debug()
patternClassifierLH = xcal.PatternClassifier([x, y//2], patternClassifierLH = xcal.PatternClassifier([x, y//2],
noiseMap[:,:y//2], noiseMap[:, :y//2],
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 = 3, tagFirstSingles=3,
nCells=memoryCells, nCells=memoryCells,
cores=cpuCores, cores=cpuCores,
allowElongated = False, allowElongated=False,
blockSize=[x, y//2], blockSize=[x, y//2],
runParallel=True) runParallel=True)
patternClassifierUH = xcal.PatternClassifier([x, y//2], patternClassifierUH = xcal.PatternClassifier([x, y//2],
noiseMap[:, y//2:], noiseMap[:, y//2:],
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 = 4, tagFirstSingles=4,
nCells=memoryCells, nCells=memoryCells,
cores=cpuCores, cores=cpuCores,
allowElongated = False, allowElongated=False,
blockSize=[x, y//2], blockSize=[x, y//2],
runParallel=True) runParallel=True)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#*****************Histogram Calculators******************# #*****************Histogram Calculators******************#
histCalOffsetCor = xcal.HistogramCalculator([x, y], histCalOffsetCor = xcal.HistogramCalculator([x, y],
bins=500, bins=500,
range=[-50, 1000], range=[-50, 1000],
nCells=memoryCells, nCells=memoryCells,
cores=cpuCores, cores=cpuCores,
blockSize=blockSize) blockSize=blockSize)
histCalPcorr = xcal.HistogramCalculator([x, y], histCalPcorr = xcal.HistogramCalculator([x, y],
bins=500, bins=500,
range=[-50, 1000], range=[-50, 1000],
nCells=memoryCells, nCells=memoryCells,
cores=cpuCores, cores=cpuCores,
blockSize=blockSize) blockSize=blockSize)
histCalPcorrS = xcal.HistogramCalculator([x, y], histCalPcorrS = xcal.HistogramCalculator([x, y],
bins=500, bins=500,
range=[-50, 1000], range=[-50, 1000],
nCells=memoryCells, nCells=memoryCells,
cores=cpuCores, cores=cpuCores,
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
patternClassifierLH._imagesPerChunk = 500 patternClassifierLH._imagesPerChunk = 500
patternClassifierUH._imagesPerChunk = 500 patternClassifierUH._imagesPerChunk = 500
#patternClassifierLH.debug() #patternClassifierLH.debug()
#patternClassifierUH.debug() #patternClassifierUH.debug()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
histCalOffsetCor.debug() histCalOffsetCor.debug()
histCalPcorr.debug() histCalPcorr.debug()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def copy_and_sanitize_non_cal_data(infile, outfile, h5base): def copy_and_sanitize_non_cal_data(infile, outfile, h5base):
if h5base.startswith("/"): if h5base.startswith("/"):
h5base = h5base[1:] h5base = h5base[1:]
dont_copy = ['image'] dont_copy = ['image']
dont_copy = [h5base+"/{}".format(do) dont_copy = [h5base+"/{}".format(do)
for do in dont_copy] for do in dont_copy]
def visitor(k, item): def visitor(k, item):
if k not in dont_copy: if k not in dont_copy:
if isinstance(item, h5py.Group): if isinstance(item, h5py.Group):
outfile.create_group(k) outfile.create_group(k)
elif isinstance(item, h5py.Dataset): elif isinstance(item, h5py.Dataset):
group = str(k).split("/") group = str(k).split("/")
group = "/".join(group[:-1]) group = "/".join(group[:-1])
infile.copy(k, outfile[group]) infile.copy(k, outfile[group])
infile.visititems(visitor) infile.visititems(visitor)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
mean_im = None mean_im = None
single_im = None single_im = None
mean_im_cc = None mean_im_cc = None
single_im_cc = None single_im_cc = None
drift_lh = [] drift_lh = []
drift_uh = [] drift_uh = []
offsetMap = np.squeeze(offsetMap) offsetMap = np.squeeze(offsetMap)
noiseMap = np.squeeze(noiseMap) noiseMap = np.squeeze(noiseMap)
badPixelMap = np.squeeze(badPixelMap) badPixelMap = np.squeeze(badPixelMap)
relGain = np.squeeze(relGain) relGain = np.squeeze(relGain)
for k, f in enumerate(file_list): for k, f in enumerate(file_list):
with h5py.File(f, 'r', driver='core') as infile: with h5py.File(f, 'r', driver='core') as infile:
out_fileb = "{}/{}".format(out_folder, f.split("/")[-1]) out_fileb = "{}/{}".format(out_folder, f.split("/")[-1])
out_file = out_fileb.replace("RAW", "CORR") out_file = out_fileb.replace("RAW", "CORR")
#out_filed = out_fileb.replace("RAW", "CORR-SC") #out_filed = out_fileb.replace("RAW", "CORR-SC")
data = None data = None
noise = None noise = None
try: try:
with h5py.File(out_file, "w") as ofile: with h5py.File(out_file, "w") as ofile:
copy_and_sanitize_non_cal_data(infile, ofile, h5path) copy_and_sanitize_non_cal_data(infile, ofile, h5path)
data = infile[h5path+"/image"][()] data = infile[h5path+"/image"][()]
nzidx = np.count_nonzero(data, axis=(1,2)) nzidx = np.count_nonzero(data, axis=(1,2))
data = data[nzidx != 0, ...] data = data[nzidx != 0, ...]
if limit_images > 0: if limit_images > 0:
data = data[:limit_images,...] data = data[:limit_images,...]
oshape = data.shape oshape = data.shape
data = np.moveaxis(data, 0, 2) data = np.moveaxis(data, 0, 2)
ddset = ofile.create_dataset(h5path+"/pixels", ddset = ofile.create_dataset(h5path+"/pixels",
oshape, oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]), chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.float32) dtype=np.float32)
ddsetm = ofile.create_dataset(h5path+"/mask", ddsetm = ofile.create_dataset(h5path+"/mask",
oshape, oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]), chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.uint32, compression="gzip") dtype=np.uint32, compression="gzip")
data = data.astype(np.float32) data = data.astype(np.float32)
offset = np.repeat(offsetMap[...,None], data.shape[2], axis=2) offset = np.repeat(offsetMap[...,None], data.shape[2], axis=2)
rg = np.repeat(relGain[:,:,None], data.shape[2], axis=2) rg = np.repeat(relGain[:,:,None], data.shape[2], axis=2)
noise = np.repeat(noiseMap[...,None], data.shape[2], axis=2) noise = np.repeat(noiseMap[...,None], data.shape[2], axis=2)
bpix = np.repeat(badPixelMap[...,None], data.shape[2], axis=2) bpix = np.repeat(badPixelMap[...,None], data.shape[2], axis=2)
data -= offset data -= offset
data *= rg data *= rg
histCalOffsetCor.fill(data) histCalOffsetCor.fill(data)
ddset[...] = np.moveaxis(data, 2, 0) ddset[...] = np.moveaxis(data, 2, 0)
ddsetm[...] = np.moveaxis(bpix, 2, 0) ddsetm[...] = np.moveaxis(bpix, 2, 0)
ofile.flush() ofile.flush()
if mean_im is None: if mean_im is None:
mean_im = np.nanmean(data, axis=2) mean_im = np.nanmean(data, axis=2)
single_im = data[...,0] single_im = data[...,0]
if do_pattern_classification: if do_pattern_classification:
#with h5py.File(out_file, "a") as ofiled: #with h5py.File(out_file, "a") as ofiled:
#copy_and_sanitize_non_cal_data(infile, ofiled, h5path) #copy_and_sanitize_non_cal_data(infile, ofiled, h5path)
ddsetcm = ofile.create_dataset(h5path+"/pixels_cm", ddsetcm = ofile.create_dataset(h5path+"/pixels_cm",
oshape, oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]), chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.float32) dtype=np.float32)
ddsetc = ofile.create_dataset(h5path+"/pixels_classified", ddsetc = ofile.create_dataset(h5path+"/pixels_classified",
oshape, oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]), chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.float32, compression="gzip") dtype=np.float32, compression="gzip")
ddsetp = ofile.create_dataset(h5path+"/patterns", ddsetp = ofile.create_dataset(h5path+"/patterns",
oshape, oshape,
chunks=(chunk_size_idim, oshape[1], oshape[2]), chunks=(chunk_size_idim, oshape[1], oshape[2]),
dtype=np.int32, compression="gzip") dtype=np.int32, compression="gzip")
patternClassifierLH._noisemap = noise[:, :x//2, :] patternClassifierLH._noisemap = noise[:, :x//2, :]
patternClassifierUH._noisemap = noise[:, x//2:, :] patternClassifierUH._noisemap = noise[:, x//2:, :]
data = cmCorrection.correct(data) data = cmCorrection.correct(data.astype(np.float32),
cellTable=np.zeros(data.shape[2], np.int32))
ddsetcm[...] = np.moveaxis(data, 2, 0) ddsetcm[...] = np.moveaxis(data, 2, 0)
dataLH = data[:, :x//2, :] dataLH = data[:, :x//2, :]
dataUH = data[:, x//2:, :] dataUH = data[:, x//2:, :]
dataLH, patternsLH = patternClassifierLH.classify(dataLH) dataLH, patternsLH = patternClassifierLH.classify(dataLH)
dataUH, patternsUH = patternClassifierUH.classify(dataUH) dataUH, patternsUH = patternClassifierUH.classify(dataUH)
data[:, :x//2, :] = dataLH data[:, :x//2, :] = dataLH
data[:, x//2:, :] = dataUH data[:, x//2:, :] = dataUH
patterns = np.zeros(data.shape, patternsLH.dtype) patterns = np.zeros(data.shape, patternsLH.dtype)
patterns[:, :x//2, :] = patternsLH patterns[:, :x//2, :] = patternsLH
patterns[:, x//2:, :] = patternsUH patterns[:, x//2:, :] = patternsUH
data[data < split_evt_primary_threshold*noise] = 0 data[data < split_evt_primary_threshold*noise] = 0
ddsetc[...] = np.moveaxis(data, 2, 0) ddsetc[...] = np.moveaxis(data, 2, 0)
ddsetp[...] = np.moveaxis(patterns, 2, 0) ddsetp[...] = np.moveaxis(patterns, 2, 0)
histCalPcorr.fill(data) histCalPcorr.fill(data)
data[patterns != 100] = np.nan data[patterns != 100] = np.nan
histCalPcorrS.fill(data) histCalPcorrS.fill(data)
if mean_im_cc is None: if mean_im_cc is None:
mean_im_cc = np.nanmean(data, axis=2) mean_im_cc = np.nanmean(data, axis=2)
single_im_cc = data[...,0] single_im_cc = data[...,0]
except Exception as e: except Exception as e:
print("Couldn't calibrate data in {}: {}".format(f, e)) print("Couldn't calibrate data in {}: {}".format(f, e))
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Common mode statistics ## ## Common mode statistics ##
Common mode statistics of at most the first 500 frames. Frames increase upwards on the y-axis, columns in the image are on the x-axis. The four quadrants are shown individually- Common mode statistics of at most the first 500 frames. Frames increase upwards on the y-axis, columns in the image are on the x-axis. The four quadrants are shown individually-
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
cmStats = cmCorrection.get() cmStats = cmCorrection.get()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
for i in range(0, cmStats.shape[0]): for i in range(0, cmStats.shape[0]):
cmStatConcat = cmStats[i, :, :500] cmStatConcat = cmStats[i, :, :500]
fig = xana.heatmapPlot(cmStatConcat.T, figsize=(8,14), lut_label="Common mode (ADU)", fig = xana.heatmapPlot(cmStatConcat.T, figsize=(8,14), lut_label="Common mode (ADU)",
y_label = "frame", y_label = "frame",
x_label = "column", x_label = "column",
vmin=-5000, vmax=5000) vmin=-5000, vmax=5000)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if do_pattern_classification: if do_pattern_classification:
print("******************LOWER HEMISPHERE******************\n") print("******************LOWER HEMISPHERE******************\n")
patternStatsLH = patternClassifierLH.getPatternStats() patternStatsLH = patternClassifierLH.getPatternStats()
fig = plt.figure(figsize=(15,15)) fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(4,4,1) ax = fig.add_subplot(4,4,1)
sfields = ["singles", "first singles", "clusters"] sfields = ["singles", "first singles", "clusters"]
mfields = ["doubles", "triples", "quads"] mfields = ["doubles", "triples", "quads"]
relativeOccurances = [] relativeOccurances = []
labels = [] labels = []
for i, f in enumerate(sfields): for i, f in enumerate(sfields):
relativeOccurances.append(patternStatsLH[f]) relativeOccurances.append(patternStatsLH[f])
labels.append(f) labels.append(f)
for i, f in enumerate(mfields): for i, f in enumerate(mfields):
for k in range(len(patternStatsLH[f])): for k in range(len(patternStatsLH[f])):
relativeOccurances.append(patternStatsLH[f][k]) relativeOccurances.append(patternStatsLH[f][k])
labels.append(f+"("+str(k)+")") labels.append(f+"("+str(k)+")")
relativeOccurances = np.array(relativeOccurances, np.float) relativeOccurances = np.array(relativeOccurances, np.float)
relativeOccurances/=np.sum(relativeOccurances) relativeOccurances/=np.sum(relativeOccurances)
pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True) pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)
ax.set_title("Pattern occurrence") ax.set_title("Pattern occurrence")
# Set aspect ratio to be equal so that pie is drawn as a circle. # Set aspect ratio to be equal so that pie is drawn as a circle.
a = ax.axis('equal') a = ax.axis('equal')
smaps = ["singlemap", "firstsinglemap", "clustermap"] smaps = ["singlemap", "firstsinglemap", "clustermap"]
for i, m in enumerate(smaps): for i, m in enumerate(smaps):
ax = fig.add_subplot(4,4,2+i) ax = fig.add_subplot(4,4,2+i)
pmap = ax.imshow(patternStatsLH[m], interpolation="nearest", vmax=2*np.nanmedian(patternStatsLH[m])) pmap = ax.imshow(patternStatsLH[m], interpolation="nearest", vmax=2*np.nanmedian(patternStatsLH[m]))
ax.set_title(m) ax.set_title(m)
cb = fig.colorbar(pmap) cb = fig.colorbar(pmap)
mmaps = ["doublemap", "triplemap", "quadmap"] mmaps = ["doublemap", "triplemap", "quadmap"]
k = 0 k = 0
for i, m in enumerate(mmaps): for i, m in enumerate(mmaps):
for j in range(4): for j in range(4):
ax = fig.add_subplot(4,4,2+len(smaps)+k) 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])) pmap = ax.imshow(patternStatsLH[m][j], interpolation="nearest", vmax=2*np.median(patternStatsLH[m][j]))
ax.set_title(m+"("+str(j)+")") ax.set_title(m+"("+str(j)+")")
cb = fig.colorbar(pmap) cb = fig.colorbar(pmap)
k+=1 k+=1
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if do_pattern_classification: if do_pattern_classification:
patternStatsUH = patternClassifierUH.getPatternStats() patternStatsUH = patternClassifierUH.getPatternStats()
fig = plt.figure(figsize=(15,15)) fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(4,4,1) ax = fig.add_subplot(4,4,1)
sfields = ["singles", "first singles", "clusters"] sfields = ["singles", "first singles", "clusters"]
mfields = ["doubles", "triples", "quads"] mfields = ["doubles", "triples", "quads"]
relativeOccurances = [] relativeOccurances = []
labels = [] labels = []
for i, f in enumerate(sfields): for i, f in enumerate(sfields):
relativeOccurances.append(patternStatsUH[f]) relativeOccurances.append(patternStatsUH[f])
labels.append(f) labels.append(f)
for i, f in enumerate(mfields): for i, f in enumerate(mfields):
for k in range(len(patternStatsUH[f])): for k in range(len(patternStatsUH[f])):
relativeOccurances.append(patternStatsUH[f][k]) relativeOccurances.append(patternStatsUH[f][k])
labels.append(f+"("+str(k)+")") labels.append(f+"("+str(k)+")")
relativeOccurances = np.array(relativeOccurances, np.float) relativeOccurances = np.array(relativeOccurances, np.float)
relativeOccurances/=np.sum(relativeOccurances) relativeOccurances/=np.sum(relativeOccurances)
pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True) pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)
ax.set_title("Pattern occurrence") ax.set_title("Pattern occurrence")
# Set aspect ratio to be equal so that pie is drawn as a circle. # Set aspect ratio to be equal so that pie is drawn as a circle.
a = ax.axis('equal') a = ax.axis('equal')
smaps = ["singlemap", "firstsinglemap", "clustermap"] smaps = ["singlemap", "firstsinglemap", "clustermap"]
for i, m in enumerate(smaps): for i, m in enumerate(smaps):
ax = fig.add_subplot(4,4,2+i) ax = fig.add_subplot(4,4,2+i)
pmap = ax.imshow(patternStatsUH[m], interpolation="nearest", vmax=2*np.nanmedian(patternStatsUH[m])) pmap = ax.imshow(patternStatsUH[m], interpolation="nearest", vmax=2*np.nanmedian(patternStatsUH[m]))
ax.set_title(m) ax.set_title(m)
cb = fig.colorbar(pmap) cb = fig.colorbar(pmap)
mmaps = ["doublemap", "triplemap", "quadmap"] mmaps = ["doublemap", "triplemap", "quadmap"]
k = 0 k = 0
for i, m in enumerate(mmaps): for i, m in enumerate(mmaps):
for j in range(4): for j in range(4):
ax = fig.add_subplot(4,4,2+len(smaps)+k) ax = fig.add_subplot(4,4,2+len(smaps)+k)
pmap = ax.imshow(patternStatsUH[m][j], interpolation="nearest", vmax=np.median(patternStatsUH[m][j])) pmap = ax.imshow(patternStatsUH[m][j], interpolation="nearest", vmax=np.median(patternStatsUH[m][j]))
ax.set_title(m+"("+str(j)+")") ax.set_title(m+"("+str(j)+")")
cb = fig.colorbar(pmap) cb = fig.colorbar(pmap)
k+=1 k+=1
print("******************UPPER HEMISPHERE******************\n") print("******************UPPER HEMISPHERE******************\n")
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if do_pattern_classification: if do_pattern_classification:
t0 = PrettyTable() t0 = PrettyTable()
t0.title = "Total number of Counts after all corrections" t0.title = "Total number of Counts after all corrections"
t0.field_names = ["Hemisphere","Singles", "First-singles", "Clusters"] t0.field_names = ["Hemisphere","Singles", "First-singles", "Clusters"]
t0.add_row(["LH", patternStatsLH['singles'], patternStatsLH['first singles'], patternStatsLH['clusters']]) t0.add_row(["LH", patternStatsLH['singles'], patternStatsLH['first singles'], patternStatsLH['clusters']])
t0.add_row(["UH", patternStatsUH['singles'], patternStatsUH['first singles'], patternStatsUH['clusters']]) t0.add_row(["UH", patternStatsUH['singles'], patternStatsUH['first singles'], patternStatsUH['clusters']])
print(t0) print(t0)
t1 = PrettyTable() t1 = PrettyTable()
t1.field_names = ["Index","D-LH", "D-UH", "T-LH", "T-UH", "Q-LH", "Q-UH"] 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([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([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([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]]) t1.add_row([3, patternStatsLH['doubles'][3], patternStatsUH['doubles'][3], patternStatsLH['triples'][3], patternStatsUH['triples'][3], patternStatsLH['quads'][3], patternStatsUH['quads'][3]])
print(t1) print(t1)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if do_pattern_classification: if do_pattern_classification:
doublesLH = patternStatsLH['doubles'][0] + patternStatsLH['doubles'][1] + patternStatsLH['doubles'][2] + patternStatsLH['doubles'][3] 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] 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] quadsLH = patternStatsLH['quads'][0] + patternStatsLH['quads'][1] + patternStatsLH['quads'][2] + patternStatsLH['quads'][3]
allsinglesLH = patternStatsLH['singles'] + patternStatsLH['first singles'] allsinglesLH = patternStatsLH['singles'] + patternStatsLH['first singles']
eventsLH = allsinglesLH + doublesLH + triplesLH + quadsLH eventsLH = allsinglesLH + doublesLH + triplesLH + quadsLH
doublesUH = patternStatsUH['doubles'][0] + patternStatsUH['doubles'][1] + patternStatsUH['doubles'][2] + patternStatsUH['doubles'][3] 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] 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] quadsUH = patternStatsUH['quads'][0] + patternStatsUH['quads'][1] + patternStatsUH['quads'][2] + patternStatsUH['quads'][3]
allsinglesUH = patternStatsUH['singles'] + patternStatsUH['first singles'] allsinglesUH = patternStatsUH['singles'] + patternStatsUH['first singles']
eventsUH = allsinglesUH + doublesUH + triplesUH + quadsUH eventsUH = allsinglesUH + doublesUH + triplesUH + quadsUH
if eventsLH > 0.: if eventsLH > 0.:
reloccurLH = np.array([allsinglesLH/eventsLH, doublesLH/eventsLH, triplesLH/eventsLH, quadsLH/eventsLH]) reloccurLH = np.array([allsinglesLH/eventsLH, doublesLH/eventsLH, triplesLH/eventsLH, quadsLH/eventsLH])
else: else:
reloccurLH = np.array([0]*4) reloccurLH = np.array([0]*4)
if eventsUH > 0.: if eventsUH > 0.:
reloccurUH = np.array([allsinglesUH/eventsUH, doublesUH/eventsUH, triplesUH/eventsUH, quadsUH/eventsUH]) reloccurUH = np.array([allsinglesUH/eventsUH, doublesUH/eventsUH, triplesUH/eventsUH, quadsUH/eventsUH])
else: else:
reloccurUH = np.array([0]*4) reloccurUH = np.array([0]*4)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if do_pattern_classification: if do_pattern_classification:
fig = plt.figure(figsize=(10,5)) fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(1,2,1) ax = fig.add_subplot(1,2,1)
labels = ['singles', 'doubles', 'triples', 'quads'] labels = ['singles', 'doubles', 'triples', 'quads']
pie = ax.pie(reloccurLH, labels=labels, autopct='%1.1f%%', shadow=True) pie = ax.pie(reloccurLH, labels=labels, autopct='%1.1f%%', shadow=True)
ax.set_title("Pattern occurrence LH") ax.set_title("Pattern occurrence LH")
# Set aspect ratio to be equal so that pie is drawn as a circle. # Set aspect ratio to be equal so that pie is drawn as a circle.
a = ax.axis('equal') a = ax.axis('equal')
ax = fig.add_subplot(1,2,2) ax = fig.add_subplot(1,2,2)
pie = ax.pie(reloccurUH, labels=labels, autopct='%1.1f%%', shadow=True) pie = ax.pie(reloccurUH, labels=labels, autopct='%1.1f%%', shadow=True)
ax.set_title("Pattern occurrence UH") ax.set_title("Pattern occurrence UH")
# Set aspect ratio to be equal so that pie is drawn as a circle. # Set aspect ratio to be equal so that pie is drawn as a circle.
a = ax.axis('equal') a = ax.axis('equal')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Mean Image of first Sequence ## ## Mean Image of first Sequence ##
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = xana.heatmapPlot(mean_im, fig = xana.heatmapPlot(mean_im,
x_label='Columns', y_label='Rows', x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)', lut_label='Signal (ADU)',
x_range=(0,y), x_range=(0,y),
y_range=(0,x), vmin=-50, vmax=70000) y_range=(0,x), vmin=-50, vmax=70000)
if do_pattern_classification: if do_pattern_classification:
fig = xana.heatmapPlot(mean_im_cc, fig = xana.heatmapPlot(mean_im_cc,
x_label='Columns', y_label='Rows', x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)', lut_label='Signal (ADU)',
x_range=(0,y), x_range=(0,y),
y_range=(0,x), vmin=-50, vmax=70000) y_range=(0,x), vmin=-50, vmax=70000)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Single Shot of first Sequnce ## ## Single Shot of first Sequnce ##
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = xana.heatmapPlot(single_im, fig = xana.heatmapPlot(single_im,
x_label='Columns', y_label='Rows', x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)', lut_label='Signal (ADU)',
x_range=(0,y), x_range=(0,y),
y_range=(0,x), vmin=-50, vmax=70000) y_range=(0,x), vmin=-50, vmax=70000)
if do_pattern_classification: if do_pattern_classification:
fig = xana.heatmapPlot(single_im_cc, fig = xana.heatmapPlot(single_im_cc,
x_label='Columns', y_label='Rows', x_label='Columns', y_label='Rows',
lut_label='Signal (ADU)', lut_label='Signal (ADU)',
x_range=(0,y), x_range=(0,y),
y_range=(0,x), vmin=-50, vmax=70000) y_range=(0,x), vmin=-50, vmax=70000)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
single = None single = None
doubles = [] doubles = []
triples = [] triples = []
quads = [] quads = []
import copy import copy
with h5py.File("{}/CORR-R{:04d}-PNCCD01-S{:05d}.h5".format(out_folder, run, sequences[0]), 'r', driver='core') as infile: with h5py.File("{}/CORR-R{:04d}-PNCCD01-S{:05d}.h5".format(out_folder, run, sequences[0]), 'r', driver='core') as infile:
data = infile[h5path+"/pixels_classified"][()].astype(np.float) data = infile[h5path+"/pixels_classified"][()].astype(np.float)
patterns = infile[h5path+"/patterns"][()].astype(np.float) patterns = infile[h5path+"/patterns"][()].astype(np.float)
i = 4 i = 4
single = copy.copy(data[...]) single = copy.copy(data[...])
single[patterns != 100] = np.nan single[patterns != 100] = np.nan
for d in range(200,204): for d in range(200,204):
double = copy.copy(data[...]) double = copy.copy(data[...])
double[patterns != d] = np.nan double[patterns != d] = np.nan
doubles.append(double) doubles.append(double)
for d in range(300,304): for d in range(300,304):
triple = copy.copy(data[...]) triple = copy.copy(data[...])
triple[patterns != d] = np.nan triple[patterns != d] = np.nan
triples.append(triple) triples.append(triple)
for d in range(400,404): for d in range(400,404):
quad = copy.copy(data[...]) quad = copy.copy(data[...])
quad[patterns != d] = np.nan quad[patterns != d] = np.nan
quads.append(quad) quads.append(quad)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
hA = 0 hA = 0
h, e = np.histogram(single.flatten(), bins=1000, range=(-500, 50000)) h, e = np.histogram(single.flatten(), bins=1000, range=(-500, 50000))
fig = plt.figure(figsize=(10,7)) fig = plt.figure(figsize=(10,7))
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
ax.step(e[:-1], h, color='blue', label='single pixel events') ax.step(e[:-1], h, color='blue', label='single pixel events')
ax.semilogy() ax.semilogy()
ax.set_xlabel("Energy (ADU)") ax.set_xlabel("Energy (ADU)")
ax.set_ylabel("Events") ax.set_ylabel("Events")
hA += h hA += h
h = 0 h = 0
for double in doubles: for double in doubles:
hd, e = np.histogram(double.flatten(), bins=1000, range=(-500, 50000)) hd, e = np.histogram(double.flatten(), bins=1000, range=(-500, 50000))
h += hd h += hd
hA += hd hA += hd
ax.step(e[:-1], h, color='red', label='double split events') ax.step(e[:-1], h, color='red', label='double split events')
h = 0 h = 0
for triple in triples: for triple in triples:
hd, e = np.histogram(triple.flatten(), bins=1000, range=(-500, 50000)) hd, e = np.histogram(triple.flatten(), bins=1000, range=(-500, 50000))
h += hd h += hd
hA += hd hA += hd
ax.step(e[:-1], h, color='green', label='triple split events') ax.step(e[:-1], h, color='green', label='triple split events')
h = 0 h = 0
for quad in quads: for quad in quads:
hd, e = np.histogram(quad.flatten(), bins=1000, range=(-500, 50000)) hd, e = np.histogram(quad.flatten(), bins=1000, range=(-500, 50000))
h += hd h += hd
hA += hd hA += hd
ax.step(e[:-1], h, color='purple', label='quadruple split events') ax.step(e[:-1], h, color='purple', label='quadruple split events')
ax.step(e[:-1], hA, color='grey', label='all events') ax.step(e[:-1], hA, color='grey', label='all events')
l = ax.legend() l = ax.legend()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` 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