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

use OO matplotlib instead of APU

parent c4d8813b
No related branches found
No related tags found
1 merge request!548plot legends for subplots properly
Hide whitespace changes
Inline Side-by-side
notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb
+ 65
− 60
View file @ 37798857
%% Cell type:markdown id: tags:
# Gain Characterization #
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/SPB/202030/p900138/scratch/karnem/r0203_r0204_v01/" # the folder to read histograms from, required
out_folder = "/gpfs/exfel/exp/SPB/202030/p900138/scratch/karnem/r0203_r0204_v01/" # the folder to output to, required
in_folder = "/gpfs/exfel/exp/SPB/202030/p900138/scratch/karnem/r0203_r0204_v01/" # the folder to read histograms from, required
out_folder = "" # the folder to output to, required
hist_file_template = "hists_m{:02d}_sum.h5" # the template to use to access histograms
modules = [10] # modules to correct, set to -1 for all, range allowed
raw_folder = "/gpfs/exfel/exp/MID/202030/p900137/raw" # Path to raw image data used to create histograms
proc_folder = "" # Path to corrected image data used to create histograms
run = 449 # of the run of image data used to create histograms
karabo_id = "MID_DET_AGIPD1M-1" # karabo karabo_id
karabo_da = ['-1'] # a list of data aggregators names, Default [-1] for selecting all data aggregators
receiver_id = "{}CH0" # inset for receiver devices
path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data
h5path = 'INSTRUMENT/{}/DET/{}:xtdf/' # path in the HDF5 file to images
h5path_idx = 'INDEX/{}/DET/{}:xtdf/' # path in the HDF5 file to images
h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP' # path to control information
karabo_id_control = "MID_IRU_AGIPD1M1" # karabo-id for control device
karabo_da_control = 'AGIPD1MCTRL00' # karabo DA for control infromation
use_dir_creation_date = True # use the creation data of the input dir for database queries
cal_db_interface = "tcp://max-exfl016:8015#8045" # the database interface to use
cal_db_timeout = 30000 # in milli seconds
local_output = True # output constants locally
db_output = False # output constants to database
# Fit parameters
peak_range = [-30, 30, 35, 70, 95, 135, 145, 220] # where to look for the peaks, [a0, b0, a1, b1, ...] exactly 8 elements
peak_width_range = [0, 30, 0, 35, 0, 40, 0, 45] # fit limits on the peak widths, [a0, b0, a1, b1, ...] exactly 8 elements
peak_norm_range = [0.0, -1, 0, -1, 0, -1, 0, -1] #
# Bad-pixel thresholds (gain evaluation error). Contribute to BadPixel bit "Gain_Evaluation_Error"
peak_lim = [-30, 30] # Limit of position of noise peak
d0_lim = [10, 80] # hard limits for distance between noise and first peak
peak_width_lim = [0.9, 1.55, 0.95, 1.65] # hard limits on the peak widths for first and second peak, in units of the noise peak. 4 parameters.
chi2_lim = [0, 3.0] # Hard limit on chi2/nDOF value
intensity_lim = 15 # Threshold on standard deviation of a histogram in ADU. Contribute to BadPixel bit "No_Entry"
gain_lim = [0.8, 1.2] # Threshold on gain in relative number. Contribute to BadPixel bit "Gain_deviation"
cell_range = [1, 3] # range of cell to be considered, [0,0] for all
pixel_range = [0, 0, 32, 32] # range of pixels x1,y1,x2,y2 to consider [0,0,512,128] for all
max_bins = 0 # Maximum number of bins to consider, 0 for all bins
batch_size = [1, 8, 8] # batch size: [cell,x,y]
fit_range = [0, 0] # range of a histogram considered for fitting in ADU. Dynamically evaluated in case [0,0]
n_peaks_fit = 4 # Number of gaussian peaks to fit including noise peak
fix_peaks = False # Fix distance between photon peaks
do_minos = False # This is additional feature of minuit to evaluate errors.
sigma_limit = 0. # If >0, repeat fit keeping only bins within mu +- sigma_limit*sigma
# Detector conditions
max_cells = 0 # number of memory cells used, set to 0 to automatically infer
bias_voltage = 300 # Bias voltage
acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine
gain_setting = 0.1 # the gain setting, use 0.1 to try to auto-determine
photon_energy = 8.05 # photon energy in keV
```
%% Cell type:code id: tags:
``` python
import glob
import os
import traceback
import warnings
from multiprocessing import Pool
import h5py
import matplotlib.pyplot as plt
import numpy as np
import sharedmem
import XFELDetAna.xfelpyanatools as xana
from cal_tools.agipdlib import get_bias_voltage
from cal_tools.agipdutils_ff import (
any_in,
fit_n_peaks,
gaussian,
gaussian_sum,
get_mask,
get_starting_parameters,
set_par_limits,
)
from cal_tools.ana_tools import get_range, save_dict_to_hdf5
from cal_tools.enums import BadPixelsFF
from iminuit import Minuit
from XFELDetAna.plotting.heatmap import heatmapPlot
from XFELDetAna.plotting.simpleplot import simplePlot
# %load_ext autotime
%matplotlib inline
warnings.filterwarnings('ignore')
```
%% Cell type:code id: tags:
``` python
peak_range = np.reshape(peak_range,(4,2))
peak_width_range = np.reshape(peak_width_range,(4,2))
peak_width_lim = np.reshape(peak_width_lim,(2,2))
peak_norm_range = [None if x == -1 else x for x in peak_norm_range]
peak_norm_range = np.reshape(peak_norm_range,(4,2))
module = modules[0]
```
%% Cell type:code id: tags:
``` python
# This is never used in this notebook and should be removed
# if bias_voltage == 0:
# # Read the bias voltage from files, if recorded.
# # If not available, make use of the historical voltage the detector is running at
# control_filename = f'{raw_folder}/r{run:04d}/RAW-R{run:04d}-{karabo_da_control}-S00000.h5'
# bias_voltage = get_bias_voltage(control_filename, karabo_id_control)
# bias_voltage = bias_voltage if bias_voltage is not None else 300
# print(f"Bias voltage: {bias_voltage}V")
```
%% Cell type:code id: tags:
``` python
def idx_gen(batch_start, batch_size):
"""
This generator iterate across pixels and memory cells starting
from batch_start until batch_start+batch_size
"""
for c_idx in range(batch_start[0], batch_start[0]+batch_size[0]):
for x_idx in range(batch_start[1], batch_start[1]+batch_size[1]):
for y_idx in range(batch_start[2], batch_start[2]+batch_size[2]):
yield(c_idx, x_idx, y_idx)
```
%% Cell type:code id: tags:
``` python
n_pixels_x = pixel_range[2]-pixel_range[0]
n_pixels_y = pixel_range[3]-pixel_range[1]
hist_data = {}
with h5py.File(f"{in_folder}/{hist_file_template.format(module)}", 'r') as hf:
hist_data['cellId'] = np.array(hf['cellId'][()])
hist_data['hRange'] = np.array(hf['hRange'][()])
hist_data['nBins'] = np.array(hf['nBins'][()])
if cell_range == [0,0]:
cell_range[1] = hist_data['cellId'].shape[0]
if max_bins == 0:
max_bins = hist_data['nBins']
hist_data['cellId'] = hist_data['cellId'][cell_range[0]:cell_range[1]]
hist_data['hist'] = np.array(hf['hist'][cell_range[0]:cell_range[1], :max_bins, :])
n_cells = cell_range[1]-cell_range[0]
hist_data['hist'] = hist_data['hist'].reshape(n_cells, max_bins, 512, 128)
hist_data['hist'] = hist_data['hist'][:,:,pixel_range[0]:pixel_range[2],pixel_range[1]:pixel_range[3]]
print(f'Data shape {hist_data["hist"].shape}')
bin_edges = np.linspace(hist_data['hRange'][0], hist_data['hRange'][1], int(hist_data['nBins']+1))
x = (bin_edges[1:] + bin_edges[:-1])[:max_bins] * 0.5
batches = []
for c_idx in range(0, n_cells, batch_size[0]):
for x_idx in range(0, n_pixels_x, batch_size[1]):
for y_idx in range(0, n_pixels_y, batch_size[2]):
batches.append([c_idx,x_idx,y_idx])
print(f'Number of batches {len(batches)}')
```
%% Cell type:code id: tags:
``` python
def fit_batch(batch_start):
current_result = {}
prev = None
for c_idx, x_idx, y_idx in idx_gen(batch_start, batch_size):
try:
y = hist_data['hist'][c_idx, :, x_idx, y_idx]
if prev is None:
prev, _ = get_starting_parameters(x, y, peak_range, n_peaks=n_peaks_fit)
if fit_range == [0, 0]:
frange = (prev[f'g0mean']-2*prev[f'g0sigma'],
prev[f'g{n_peaks_fit-1}mean'] + prev[f'g{n_peaks_fit-1}sigma'])
else:
frange = fit_range
set_par_limits(prev, peak_range, peak_norm_range,
peak_width_range, n_peaks_fit)
minuit = fit_n_peaks(x, y, prev, frange,
do_minos=do_minos, n_peaks=n_peaks_fit,
fix_d01=fix_peaks, sigma_limit=sigma_limit,)
ndof = np.rint(frange[1]-frange[0])-len(minuit.args) ## FIXME: this line is wrong if fix_peaks is True
current_result['chi2_ndof'] = minuit.fval/ndof
res = minuit.fitarg
if fix_peaks : ## set g2 and g3 mean correctly
for i in range(2,n_peaks_fit):
d = res[f'g1mean'] - res[f'g0mean']
res[f'g{i}mean'] = res[f'g0mean'] + d*i
current_result.update(res)
current_result.update(minuit.get_fmin())
fit_result['chi2_ndof'][c_idx, x_idx, y_idx] = current_result['chi2_ndof']
for key in res.keys():
if key in fit_result:
fit_result[key][c_idx, x_idx, y_idx] = res[key]
fit_result['mask'][c_idx, x_idx, y_idx] = get_mask(current_result,
peak_lim,
d0_lim, chi2_lim,
peak_width_lim)
except Exception as e:
fit_result['mask'][c_idx, x_idx,
y_idx] = BadPixelsFF.FIT_FAILED.value
print(c_idx, x_idx, y_idx, e, traceback.format_exc())
if fit_result['mask'][c_idx, x_idx, y_idx] == 0:
prev = res
else:
prev = None
```
%% Cell type:markdown id: tags:
# Single fit ##
Left plot shows starting parameters for fitting. Right plot shows result of the fit. Errors are evaluated with minos.
%% Cell type:code id: tags:
``` python
hist = hist_data['hist'][1,:,1, 1]
prev, shapes = get_starting_parameters(x, hist, peak_range, n_peaks=n_peaks_fit)
if fit_range == [0, 0]:
frange = (prev[f'g0mean']-2*prev[f'g0sigma'],
prev[f'g3mean'] + prev[f'g3sigma'])
else:
frange = fit_range
set_par_limits(prev, peak_range, peak_norm_range,
peak_width_range, n_peaks=n_peaks_fit)
minuit = fit_n_peaks(x, hist, prev, frange,
do_minos=True, n_peaks=n_peaks_fit,
fix_d01=fix_peaks,
sigma_limit=sigma_limit,
)
print (minuit.get_fmin())
minuit.print_matrix()
print(minuit.get_param_states())
```
%% Cell type:code id: tags:
``` python
res = minuit.fitarg
if fix_peaks :
for i in range(2,n_peaks_fit):
d = res[f'g1mean'] - res[f'g0mean']
res[f'g{i}mean'] = res[f'g0mean'] + d*i
err = minuit.errors
p = minuit.args
ya = np.arange(0,1e4)
y = gaussian_sum(x,n_peaks_fit, *p)
peak_colors = ['g', 'y', 'b', 'orange']
peak_hist = hist.copy()
d=[]
if sigma_limit > 0 :
sel2 = (np.abs(x - res['g0mean']) < sigma_limit*res['g0sigma']) | \
(np.abs(x - res['g1mean']) < sigma_limit*res['g1sigma']) | \
(np.abs(x - res['g2mean']) < sigma_limit*res['g2sigma']) | \
(np.abs(x - res['g3mean']) < sigma_limit*res['g3sigma'])
peak_hist[~sel2] = 0
valley_hist = hist.copy()
valley_hist[sel2] = 0
d.append({'x': x,
'y': valley_hist.astype(np.float64),
'y_err': np.sqrt(valley_hist),
'drawstyle': 'bars',
'errorstyle': 'bars',
'transparency': '95%',
'errorcoarsing': 3,
'label': f'X-ray Data)'
})
htitle = f'X-ray Data, (μ±{sigma_limit:0.1f}σ)'
else :
htitle = 'X-ray Data'
d.append({'x': x,
'y': peak_hist.astype(np.float64),
'y_err': np.sqrt(peak_hist),
'drawstyle': 'bars',
'errorstyle': 'bars',
'errorcoarsing': 3,
'label': htitle,
}
)
d.append({'x': x,
'y': y,
'y2': (hist-y)/np.sqrt(hist),
'drawstyle':'line',
'drawstyle2': 'steps-mid',
'label': 'Fit'
}
)
for i in range(n_peaks_fit):
d.append({'x': x,
'y': gaussian(x, res[f'g{i}n'], res[f'g{i}mean'], res[f'g{i}sigma']),
'drawstyle':'line',
'color': peak_colors[i],
})
d.append({'x': np.full_like(ya, res[f'g{i}mean']),
'y': ya,
'drawstyle': 'line',
'linestyle': 'dashed',
'color': peak_colors[i],
'label': f'peak {i} = {res[f"g{i}mean"]:0.1f} $ \pm $ {err[f"g{i}mean"]:0.2f} ADU' })
```
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(10,5))
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_size_inches(16, 7)
for i, shape in enumerate(shapes):
idx = shape[3]
plt.errorbar(x[idx], hist[idx], np.sqrt(hist[idx]),
marker='+', ls=''
)
ax1.errorbar(
x[idx], hist[idx],
np.sqrt(hist[idx]),
marker='+', ls='',
)
yg = gaussian(x[idx], *shape[:3])
l = f'Peak {i}: {shape[1]:0.1f} $ \pm $ {shape[2]:0.2f} ADU'
plt.plot(x[idx], yg, label=l)
plt.grid(True)
plt.xlabel("Signal [ADU]")
plt.ylabel("Counts")
plt.legend(ncol=2)
ax1.plot(x[idx], yg, label=l)
ax1.grid(True)
ax1.set_xlabel("Signal [ADU]")
ax1.set_ylabel("Counts")
ax1.legend(ncol=2)
plt.show()
```
%% Cell type:code id: tags:
``` python
_ = xana.simplePlot(
d,
use_axis=ax2,
x_label='Signal [ADU]',
y_label='Counts',
secondpanel=True, y_log=False,
x_range=(frange[0], frange[1]),
y_range=(1., np.max(hist)*1.6),
legend='top-left-frame-ncol2'
legend='top-left-frame-ncol2',
)
plt.show()
```
%% Cell type:markdown id: tags:
## All fits ##
%% Cell type:code id: tags:
``` python
# Allocate memory for fit results
fit_result = {}
keys = list(minuit.fitarg.keys())
keys = [x for x in keys if 'limit_' not in x and 'fix_' not in x]
keys += ['chi2_ndof', 'mask', 'gain']
for key in keys:
dtype = 'f4'
if key == 'mask':
dtype = 'i4'
fit_result[key] = sharedmem.empty([n_cells, n_pixels_x, n_pixels_y], dtype=dtype)
```
%% Cell type:code id: tags:
``` python
# Perform fitting
with Pool() as pool:
const_out = pool.map(fit_batch, batches)
```
%% Cell type:code id: tags:
``` python
# Evaluate bad pixels
fit_result['gain'] = (fit_result['g1mean'] - fit_result['g0mean'])/photon_energy
# Calculate histogram width and evaluate cut
h_sums = np.sum(hist_data['hist'], axis=1)
hist_norm = hist_data['hist'] / h_sums[:, None, :, :]
hist_mean = np.sum(hist_norm[:, :max_bins, ...] *
x[None, :, None, None], axis=1)
hist_sqr = (x[None, :, None, None] - hist_mean[:, None, ...])**2
hist_std = np.sqrt(np.sum(hist_norm[:, :max_bins, ...] * hist_sqr, axis=1))
fit_result['mask'][hist_std<intensity_lim] |= BadPixelsFF.NO_ENTRY.value
# Bad pixel on gain deviation
gains = np.copy(fit_result['gain'])
gains[fit_result['mask']>0] = np.nan
gain_mean = np.nanmean(gains, axis=(1,2))
fit_result['mask'][fit_result['gain'] > gain_mean[:,None,None]*gain_lim[1] ] |= BadPixelsFF.GAIN_DEVIATION.value
fit_result['mask'][fit_result['gain'] < gain_mean[:,None,None]*gain_lim[0] ] |= BadPixelsFF.GAIN_DEVIATION.value
```
%% Cell type:code id: tags:
``` python
# Save fit results
os.makedirs(out_folder, exist_ok=True)
out_name = f'{out_folder}/fits_m{module:02d}.h5'
print(f'Save to file: {out_name}')
save_dict_to_hdf5({'data': fit_result}, out_name)
```
%% Cell type:markdown id: tags:
## Summary across cells ##
%% Cell type:code id: tags:
``` python
labels = ['Noise peak [ADU]',
'First photon peak [ADU]',
f"gain [ADU/keV], $\gamma$={photon_energy} [keV]",
"$\chi^2$/nDOF",
'Fraction of bad pixels']
for i, key in enumerate(['g0mean', 'g1mean', 'gain', 'chi2_ndof', 'mask']):
labels = [
"Noise peak [ADU]",
"First photon peak [ADU]",
f"gain [ADU/keV] $\gamma$={photon_energy} [keV]",
"$\chi^2$/nDOF",
"Fraction of bad pixels",
]
for i, key in enumerate(['g0mean', 'g1mean', 'gain', 'chi2_ndof', 'mask']):
fig = plt.figure(figsize=(20,5))
ax = fig.add_subplot(121)
data = fit_result[key]
if key == 'mask':
data = data>0
data = data > 0
vmin, vmax = [0, 1]
else:
vmin, vmax = get_range(data, 5)
_ = heatmapPlot(np.mean(data, axis=0).T,
add_panels=False, cmap='viridis', use_axis=ax,
vmin=vmin, vmax=vmax, lut_label=labels[i] )
_ = heatmapPlot(
np.mean(data, axis=0).T,
add_panels=False, cmap='viridis', use_axis=ax,
vmin=vmin, vmax=vmax, lut_label=labels[i]
)
if key != 'mask':
vmin, vmax = get_range(data, 7)
ax1 = fig.add_subplot(122)
_ = xana.histPlot(ax1,data.flatten(),
bins=45,range=[vmin, vmax],
log=True,color='red',histtype='stepfilled')
plt.xlabel(labels[i])
plt.ylabel("Counts")
ax = fig.add_subplot(122)
_ = xana.histPlot(
ax, data.flatten(),
bins=45,range=[vmin, vmax],
log=True,color='red',histtype='stepfilled'
)
ax.set_xlabel(labels[i])
ax.set_ylabel("Counts")
```
%% Cell type:markdown id: tags:
## histograms of fit parameters ##
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111)
a = ax.hist(hist_std.flatten(), bins=100, range=(0,100) )
ax.plot([intensity_lim, intensity_lim], [0, np.nanmax(a[0])], linewidth=1.5, color='red' )
ax.set_xlabel('Histogram width [ADU]', fontsize=14)
ax.set_ylabel('Number of histograms', fontsize=14)
ax.set_title(f'{hist_std[hist_std<intensity_lim].shape[0]} histograms below threshold in {intensity_lim} ADU',
fontsize=14, fontweight='bold')
ax.grid()
plt.yscale('log')
ax.set_yscale('log')
```
%% Cell type:code id: tags:
``` python
def plot_par_distr(par):
fig = plt.figure(figsize=(16,5))
fig = plt.figure(figsize=(16, 5))
sel = fit_result['mask'] == 0
for i in range(n_peaks_fit) :
data=fit_result[f"g{i}{par}"]
plt_range=(-1,50)
if par =='mean':
plt_range=[peak_range[i][0] ,peak_range[i][1]]
num_bins = int(plt_range[1] - plt_range[0])
ax = fig.add_subplot(1,n_peaks_fit,i+1)
_ = xana.histPlot(ax,data.flatten(),
bins= num_bins,range=plt_range,
log=True,color='red',
label='all fits',)
a = ax.hist(data[sel].flatten(),
bins=num_bins, range=plt_range,
log=True,color='g',
label='good fits only',
)
plt.xlabel(f"g{i} {par} [ADU]")
ax.set_xlabel(f"g{i} {par} [ADU]")
ax.legend()
plot_par_distr('mean')
plot_par_distr('sigma')
```
%% Cell type:code id: tags:
``` python
sel = fit_result['mask'] == 0
dsets = {'d01 [ADU]':fit_result[f"g1mean"]-fit_result[f"g0mean"],
'gain [ADU/keV]':fit_result[f"gain"],
'gain relative to module mean':fit_result[f"gain"]/np.mean(gain_mean),
}
fig = plt.figure(figsize=(16,5))
for i, (par, data) in enumerate(dsets.items()):
ax = fig.add_subplot(1,3,i+1)
ax = fig.add_subplot(1, 3, i+1)
plt_range=get_range(data, 10)
num_bins = 100
_ = xana.histPlot(ax,data.flatten(),
bins= num_bins,range=plt_range,
log=True,color='red',
label='all fits',)
a = ax.hist(data[sel].flatten(),
bins=num_bins, range=plt_range,
log=True,color='g',
label='good fits only',
)
plt.xlabel(f"{par}")
ax.set_xlabel(f"{par}")
ax.legend()
if 'd01' in par :
plt.axvline(d0_lim[0])
plt.axvline(d0_lim[1])
ax.axvline(d0_lim[0])
ax.axvline(d0_lim[1])
if 'rel' in par :
plt.axvline(gain_lim[0])
plt.axvline(gain_lim[1])
ax.axvline(gain_lim[0])
ax.axvline(gain_lim[1])
```
%% Cell type:markdown id: tags:
## Summary across pixels ##
Mean and median values are calculated across all pixels for each memory cell.
%% Cell type:code id: tags:
``` python
def plot_error_band(key, x, ax):
cdata = np.copy(fit_result[key])
cdata[fit_result['mask']>0] = np.nan
mean = np.nanmean(cdata, axis=(1,2))
median = np.nanmedian(cdata, axis=(1,2))
std = np.nanstd(cdata, axis=(1,2))
mad = np.nanmedian(np.abs(cdata - median[:,None,None]), axis=(1,2))
ax.plot(x, mean, 'k', color='#3F7F4C', label=" mean value ")
ax.plot(x, median, 'o', color='red', label=" median value ")
ax.fill_between(x, mean-std, mean+std,
alpha=0.6, edgecolor='#3F7F4C', facecolor='#7EFF99',
linewidth=1, linestyle='dashdot', antialiased=True,
label=" mean value $ \pm $ std ")
ax.fill_between(x, median-mad, median+mad,
alpha=0.3, edgecolor='red', facecolor='red',
linewidth=1, linestyle='dashdot', antialiased=True,
label=" median value $ \pm $ mad ")
if f'error_{key}' in fit_result:
cerr = np.copy(fit_result[f'error_{key}'])
cerr[fit_result['mask']>0] = np.nan
meanerr = np.nanmean(cerr, axis=(1,2))
ax.fill_between(x, mean-meanerr, mean+meanerr,
alpha=0.6, edgecolor='#089FFF', facecolor='#089FFF',
linewidth=1, linestyle='dashdot', antialiased=True,
label=" mean fit error ")
x = np.linspace(*cell_range, n_cells)
for i, key in enumerate(['g0mean', 'g1mean', 'gain', 'chi2_ndof']):
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111)
plot_error_band(key, x, ax)
ax.set_xlabel('Memory Cell ID', fontsize=14)
ax.set_ylabel(labels[i], fontsize=14)
ax.grid()
ax.legend()
```
%% Cell type:markdown id: tags:
## Cut flow ##
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(10, 5))
fig, ax = plt.subplots()
fig.set_size_inches(10, 5)
n_bars = 8
x = np.arange(n_bars)
width = 0.3
msk = fit_result['mask']
n_fits = np.prod(msk.shape)
y = [any_in(msk, BadPixelsFF.FIT_FAILED.value),
any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value),
any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
BadPixelsFF.CHI2_THRESHOLD.value),
any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value),
any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |
BadPixelsFF.NOISE_PEAK_THRESHOLD.value),
any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |
BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value),
any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |
BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value
| BadPixelsFF.NO_ENTRY.value),
any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |
BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value
| BadPixelsFF.NO_ENTRY.value| BadPixelsFF.GAIN_DEVIATION.value)
]
y2 = [any_in(msk, BadPixelsFF.FIT_FAILED.value),
any_in(msk, BadPixelsFF.ACCURATE_COVAR.value),
any_in(msk, BadPixelsFF.CHI2_THRESHOLD.value),
any_in(msk, BadPixelsFF.GAIN_THRESHOLD.value),
any_in(msk, BadPixelsFF.NOISE_PEAK_THRESHOLD.value),
any_in(msk, BadPixelsFF.PEAK_WIDTH_THRESHOLD.value),
any_in(msk, BadPixelsFF.NO_ENTRY.value),
any_in(msk, BadPixelsFF.GAIN_DEVIATION.value)
]
y = (1 - np.sum(y, axis=(1,2,3))/n_fits)*100
y2 = (1 - np.sum(y2, axis=(1,2,3))/n_fits)*100
labels = ['Fit failes',
'Accurate covar',
'Chi2/nDOF',
'Gain',
'Noise peak',
'Peak width',
'No Entry',
'Gain deviation']
plt.bar(x, y2, width, label='Only this cut')
plt.bar(x, y, width, label='Cut flow')
plt.xticks(x, labels, rotation=90)
plt.ylim(y[5]-0.5, 100)
plt.grid(True)
plt.legend()
ax.bar(x, y2, width, label='Only this cut')
ax.bar(x, y, width, label='Cut flow')
ax.set_xticks(x)
ax.set_xticklabels(labels, rotation=90)
ax.set_ylim(y[5]-0.5, 100)
ax.grid(True)
ax.legend()
plt.show()
```
......
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