Stream Analysis¶
This notebook provides an analysis of the software triggered data of a detector module. This includes the calculation of rate and stability cuts, energy spectra and light yield plots, and efficiencies.
print('Let´s start!')
Let´s start!
Introduction¶
This notebook builds upon the hardware data analysis notebook and the trigger script. We do cuts and the energy calibration.
First we import packages.
import cait as ai
import matplotlib.pyplot as plt
import numpy as np
import h5py
import pickle
from scipy.stats import norm, expon
%config InlineBackend.figure_formats = ['svg'] # we need this for a suitable resolution of the plots
And we define a set of constants and paths.
RUN = ... # put an string for the number of the experiments run, e.g. '34'
MODULE = ... # put a name for the detector, e.g. 'DetA'
PATH_PROC_DATA = ... # path to where you want to store the HDF5 files
FILE_NMBRS = [] # a list of string, the file number you want to analyse, e.g. ['001', '002', '003']
RDT_CHANNELS = [] # a list of strings of the channels, e.g. [0, 1] (written in PAR file - attention, the PAR file counts from 1, Cait from 0)
RECORD_LENGTH = 16384 # the number of samples within a record window (read in PAR file)
SAMPLE_FREQUENCY = 25000 # the sample frequency of the measurement (read in PAR file)
DOWN_SEF = 4 # the downsample rate for the standard event fit
DOWN_BLF = 16 # the downsample rate for the baseline fit
PROCESSES = 8 # the number of processes for parallelization
PCA_COMPONENTS = 2 # the number of pca components to calculate
SKIP_FNMR = [] # in case the loop crashed at some point and you want to start from a specific file number, write here the numbers to ignore, e.g. ['001', '002']
# typically you do not need to change the values below
FNAME_EFFICIENCY = 'efficiency_{:03d}'.format(len(FILE_NMBRS) - 1)
FNAME_TRAINING = 'training_001'
FNAME_STREAM = 'stream_{:03d}'.format(len(FILE_NMBRS) - 1)
H5_CHANNELS = list(range(len(RDT_CHANNELS)))
SEF_APP = '_down{}'.format(DOWN_SEF) if DOWN_SEF > 1 else ''
We need some values that we already calculated from the hardware data.
DETECTOR_MASS = ... # the detector mass in kg
print('Detector mass in kg: ', DETECTOR_MASS)
BL_RESOLUTION_OF = [] # list of the baseline resolutions, calculated with the superposition method, in mV
THRESHOLDS = [(6.5 * r) * 1e-3 for r in BL_RESOLUTION_OF]
print('OF resolution in V: ', BL_RESOLUTION_OF)
print('OF thresholds in V: ', THRESHOLDS)
FIT_BL_SIGMAS = [] # list of the baseline resolutions, calculated with the noise fit model, in V
FIT_THRESHOLDS = [] # list of the trigger thresholds, calculated with the noise fit model, in V
print('Fit thresholds in V: ', FIT_THRESHOLDS)
TRUNCATION_LEVELS = [] # list of the truncation levels
MAXIMAL_EVENT_HEIGHTS = [] # list of the maximal event heights to be included in the energy calibration
PATH_PULSER_MODEL = ... # put an string of the path where you want to store the pulser model
PEAK_ENERGY = 5.89 * 8/9 + 6.49 * 1/9 # in case of an iron source - change this for different source
Here we assemble calculated values from further down in the notebook. Fill them up while you go.
PEC_FACTOR = [...] # list of the values source_peak_energy/source_peak_pulse_height
print('Energy resolution (eV): ', [b*p for b,p in zip(BL_RESOLUTION_OF, PEC_FACTOR)])
print('Treshold (eV): ', [1000*b*p for b,p in zip(THRESHOLDS, PEC_FACTOR)])
print('Treshold (eV) 1 Noise Trigger Fit: ', [1000*b*p for b,p in zip(FIT_THRESHOLDS, PEC_FACTOR)])
print('Treshold from Trigger Efficiency (eV): ', 1000*np.array([...])) # The fitted trigger efficiency threshold
CPE_FACTOR = [...]# list of the values source_peak_energy/source_peak_equivalent_tpa_value
LIVE_TIME = ... # the total live time of the detector, calculated at the end of the file
EXPOSURE = LIVE_TIME/24 * DETECTOR_MASS
print('Live Time in h: ', LIVE_TIME)
print('Exposure in kg * days: ', EXPOSURE)
We create a DataHandler instance to do the processing.
dh_stream = ai.DataHandler(run=RUN,
module=MODULE,
channels=RDT_CHANNELS)
dh_stream.set_filepath(path_h5=PATH_PROC_DATA,
fname=FNAME_STREAM,
appendix=False)
View Events¶
Lets have a look at the stream triggered events to check, if the triggering worked properly.
ei = ai.EventInterface(module=MODULE,
run=RUN,
nmbr_channels=len(H5_CHANNELS))
ei.load_h5(path=PATH_PROC_DATA,
fname=FNAME_STREAM,
channels=RDT_CHANNELS,
appendix=False,
which_to_label=['events'])
We need to include a labels file.
# ei.create_labels_csv(path=PATH_PROC_DATA)
# ei.load_labels_csv(path=PATH_PROC_DATA)
We can also load the optimal filter transfer function to view the events the way they were triggered.
# ei.load_of(group_name_appendix='_tp')
# ei.load_sev_par(name_appendix='_down4', group_name_appendix='_tp')
# ei.load_of()
# ei.load_sev_par(name_appendix='_down4')
Okay, lets see the events.
ei.start(start_from_idx=0, print_label_list=False)
Cuts¶
To clean the data from artifacts, we want to apply some quality cuts.
In case we have not included the SEV, NPS and OF during the trigger loop, and have not calculated the corresponding features, we can do this here.
# Only run this if it wasnt run during the trigger loop!
dh_hw = ai.DataHandler(run=RUN,
module=MODULE,
channels=RDT_CHANNELS)
dh_hw.set_filepath(path_h5=PATH_PROC_DATA,
fname='hw_{:03d}'.format(len(FILE_NMBRS)-1),
appendix=False)
# Only run this if it wasnt run during the trigger loop!
dh_stream.include_sev(sev=dh_hw.get('stdevent','event'),
fitpar=dh_hw.get('stdevent','fitpar'),
mainpar=dh_hw.get('stdevent','mainpar'))
dh_stream.include_nps(nps=dh_hw.get('noise','nps'))
dh_stream.include_of(of_real=dh_hw.get('optimumfilter','optimumfilter_real'),
of_imag=dh_hw.get('optimumfilter','optimumfilter_imag'))
dh_stream.include_sev(sev=dh_hw.get('stdevent_tp','event'),
fitpar=dh_hw.get('stdevent_tp','fitpar'),
mainpar=dh_hw.get('stdevent_tp','mainpar'),
group_name_appendix='_tp')
dh_stream.include_of(of_real=dh_hw.get('optimumfilter_tp','optimumfilter_real'),
of_imag=dh_hw.get('optimumfilter_tp','optimumfilter_imag'),
group_name_appendix='_tp')
# Only run this if it wasnt run during the trigger loop!
dh_stream.apply_of(first_channel_dominant=True)
dh_stream.apply_of(type='testpulses', name_appendix_group='_tp')
dh_stream.apply_sev_fit(down=DOWN_SEF, name_appendix='_down{}'.format(DOWN_SEF), processes=PROCESSES,
truncation_level=TRUNCATION_LEVELS, verb=True, first_channel_dominant=True)
dh_stream.apply_sev_fit(type='testpulses', group_name_appendix='_tp',
down=DOWN_SEF, name_appendix='_down{}'.format(DOWN_SEF), processes=PROCESSES,
truncation_level=TRUNCATION_LEVELS, verb=True)
# do your own cuts!
clean_events = ai.cuts.LogicalCut(initial_condition=np.abs(dh_stream.get('events', 'mainpar')[0,:,8]) < 2e-6)
clean_events.add_condition(np.abs(dh_stream.get('events', 'mainpar')[1,:,8]) < 2e-6)
clean_events.add_condition(dh_stream.get('events', 'mainpar')[0,:,0] < 1)
clean_events.add_condition(dh_stream.get('events', 'mainpar')[1,:,0] < 1.5)
clean_events.add_condition(dh_stream.get('events', 'mainpar')[0,:,3] < 4500)
clean_events.add_condition(dh_stream.get('events', 'mainpar')[0,:,3] > 3900)
clean_events.add_condition(dh_stream.get('events', 'mainpar')[1,:,3] < 4500)
clean_events.add_condition(dh_stream.get('events', 'mainpar')[1,:,3] > 3900)
for c in H5_CHANNELS:
dh_stream.apply_logical_cut(cut_flag=clean_events.get_flag(),
naming='clean_events',
channel=c,
type='events',
delete_old=False)
dh_stream.apply_pca(nmbr_components=PCA_COMPONENTS,
down=DOWN_SEF,
fit_idx=clean_events.get_idx())
One cut we want to apply is a rate cut.
dh_stream.calc_rate_cut(min=0, max=30)
Now we plot the histogram of the rate.
nmbr_bins = 500
total_hours = 1000
good_intervals = dh_stream.get('metainfo', 'rate_stable')
dh_stream.show_values(group='events', key='hours', bins=nmbr_bins,
xlabel='Time (h)', ylabel=f'Counts/({total_hours/nmbr_bins*60:.1f} min)', title='Event Rate',
show=False)
for iv in good_intervals:
plt.axvspan(xmin=iv[0], xmax=iv[1], color='green', alpha=0.3, zorder=80)
# plt.xlim((400,450))
plt.show()
Do the stability cut.
cp_ph = dh_stream.get('controlpulses', 'pulse_height')
bounds = []
for c in H5_CHANNELS:
median_cp = np.median(cp_ph[c])
print('Median pulse height Channel {}: {} V'.format(c, median_cp))
bounds.append((median_cp*0.8, median_cp*1.2))
dh_stream.calc_controlpulse_stability(channel=c, significance=3, lb=bounds[c][0], ub=bounds[c][1])
Now plot the control pulse heights.
cp_h = dh_stream.get('controlpulses', 'hours')
y_lims = [None, None]
for c in H5_CHANNELS:
instable_iv = dh_stream.get('metainfo', f'controlpulse_instable_ch{c}')
plt.close()
xlims, ylims, I = ai.styles.scatter_img(x_data=cp_h, y_data=cp_ph[c], ylims=y_lims[c], height=80, width=200, alpha=0.3)
ax_extent = list(xlims)+list(ylims)
ai.styles.use_cait_style()
plt.imshow(I,
vmin=0,
vmax=1,
cmap=plt.get_cmap('viridis'),
interpolation='lanczos',
aspect='auto',
extent=ax_extent,
)
for iv in instable_iv:
plt.axvspan(iv[0], iv[1], color='red', alpha=0.2, zorder=70)
plt.grid(alpha=0.2)
plt.title('Control Pulse Stability Channel {}'.format(c))
plt.xlabel('Time (h)')
plt.ylabel('Pulse Height (V)')
plt.show()
We do the quality cuts further down.
Energy Calibration¶
We can double check the CPE factors with the iron calibration peak, which is done below. Afterwards, we do the energy fine tuning with the test pulses.
For as estimator for the pulse height, we always take the optimum filter.
x_ranges = [(0, 1.6), (0, 0.3)]
y_ranges = [(0, 200), (0, 1000)]
for c, xr, yr in zip(H5_CHANNELS, x_ranges, y_ranges):
dh_stream.show_values(group='events', key='of_ph', bins=250, idx0=c, idx2=None, range=xr, yran=yr,
xlabel='Pulse Height (V)', ylabel='Counts', title='Spectrum PH Channel {}'.format(c), show=False)
# plt.axvline(x=calibration_peak[c], color='black', linestyle='dotted', linewidth=2.5, zorder=70) # run the cell below first!
plt.show()
ph = dh_stream.get('events','of_ph')
iv = [(0.3, 0.45), (0.008, 0.03)]
calibration_peak = []
for c in H5_CHANNELS:
peak_events = ai.cuts.LogicalCut(initial_condition=ph[c]>iv[c][0])
peak_events.add_condition(ph[c]<iv[c][1])
calibration_peak.append(np.mean(ph[c, peak_events.get_idx()]))
print('Peak position with OF: ', calibration_peak[-1])
We can here already calculate the PEC rate.
pec_factor = PEAK_ENERGY / np.array(calibration_peak)
print('The PEC factors are: ', pec_factor)
Fine Tuning¶
We will now do a fine tuning of the energy calibration. For this we first have a look at the saturation of the test pulses.
tp_tpa = dh_stream.get('testpulses', 'testpulseamplitude')
tp_ph = dh_stream.get('testpulses', 'sev_fit_par{}'.format(SEF_APP))[:, :, 0]
y_lims = [(0, 4), (0, 0.5)]
for c in H5_CHANNELS:
plt.close()
xlims, ylims, I = ai.styles.scatter_img(x_data=tp_tpa, y_data=tp_ph[c], ylims=y_lims[c], height=150, width=200, alpha=0.5)
ax_extent = list(xlims)+list(ylims)
ai.styles.use_cait_style()
plt.imshow(I,
vmin=0,
vmax=1,
cmap=plt.get_cmap('viridis'),
interpolation='lanczos',
aspect='auto',
extent=ax_extent,
)
plt.grid(alpha=0.2)
plt.title('TP Saturation Truncated Fit Channel {}'.format(c))
plt.xlabel('TPA (V)')
plt.ylabel('Pulse Height (V)')
plt.colorbar()
plt.show()
We do a stability cut on the test pulses for which we define lower and upper bounds for the pulse heights of each TPA value.
tp_tpa = dh_stream.get('testpulses', 'testpulseamplitude')
tp_ph = dh_stream.get('testpulses', 'of_ph')[:, :]
unique_tp = np.unique(tp_tpa)
print('Unique testpulse heights: ', unique_tp)
Check the plot, if the bounds look good.
idx = 1
c = 0
print('tps', unique_tp[idx])
mean = np.mean(tp_ph[c][tp_tpa == unique_tp[idx]])
std = np.std(tp_ph[c][tp_tpa == unique_tp[idx]])
med = np.median(tp_ph[c][tp_tpa == unique_tp[idx]])
quantile_low = np.quantile(tp_ph[c][tp_tpa == unique_tp[idx]], 0.18)
quantile_up = np.quantile(tp_ph[c][tp_tpa == unique_tp[idx]], 0.82)
maxi = np.max(tp_ph[c][tp_tpa == unique_tp[idx]])
mini = np.min(tp_ph[c][tp_tpa == unique_tp[idx]])
print('mean', mean)
print('std', std)
print('med', med)
print('quantile_low', quantile_low)
print('quantile_up', quantile_up)
print('maxi', maxi)
print('mini', mini)
mean_dev = quantile_up - quantile_low
plt.close()
ai.styles.use_cait_style()
plt.hist(tp_ph[c][tp_tpa == unique_tp[idx]], bins=100, range=(med - 5*mean_dev, med + 5*mean_dev), zorder=30)
plt.axvline(x=med, color='red', label='med', zorder=50)
plt.axvspan(quantile_low - 2*mean_dev, quantile_up + 2*mean_dev, alpha=0.5, color='red', zorder=15)
plt.xlabel('Pulse Height (V)')
plt.ylabel('Counts')
plt.title(f'TPA = {unique_tp[idx]}')
plt.legend()
plt.yscale('log')
ai.styles.make_grid()
plt.show()
lb = []
ub = []
for c in H5_CHANNELS:
medians = [np.median(tp_ph[c][tp_tpa == tpa]) for tpa in unique_tp]
lower_quantiles = [np.quantile(tp_ph[c][tp_tpa == tpa], 0.18) for tpa in unique_tp]
upper_quantiles = [np.quantile(tp_ph[c][tp_tpa == tpa], 0.82) for tpa in unique_tp]
mean_deviations = [u - l for l,u in zip(lower_quantiles, upper_quantiles)]
lb.append([ l - 2*m for l,m in zip(lower_quantiles, mean_deviations)])
ub.append([ u + 2*m for u,m in zip(upper_quantiles, mean_deviations)])
Check again the plot.
y_lims = [(0, 1.7), (0, 0.35)]
for c in H5_CHANNELS:
plt.close()
xlims, ylims, I = ai.styles.scatter_img(x_data=tp_tpa, y_data=tp_ph[c], ylims=y_lims[c], height=150, width=200, alpha=0.5)
ax_extent = list(xlims)+list(ylims)
ai.styles.use_cait_style()
plt.imshow(I,
vmin=0,
vmax=1,
cmap=plt.get_cmap('viridis'),
interpolation='lanczos',
aspect='auto',
extent=ax_extent,
)
plt.grid(alpha=0.2)
plt.plot(unique_tp, lb[c], color='red', linewidth=1.5)
plt.plot(unique_tp, ub[c], color='red', linewidth=1.5)
plt.title('TP Saturation OF Channel {}'.format(c))
plt.xlabel('TPA (V)')
plt.ylabel('Pulse Height (V)')
plt.colorbar()
plt.show()
Now do the test pulse stability cut.
for c in H5_CHANNELS:
dh_stream.calc_testpulse_stability(c, significance=3, ub = ub[c], lb = lb[c])
Now we are ready to do the actual energy calibration. For this we use a gradient boosted decision tree.
pm = dh_stream.calc_calibration(starts_saturation=[1.6, 0.3], # stop energy calibration at these values
cpe_factor=[1, 1],
plot=False,
only_stable=True,
exclude_tpas=[],
interpolation_method='tree',
poly_order=5,
method='of',
return_pulser_models=True,
#name_appendix_ev='_down{}'.format(down),
#name_appendix_tp='_down{}'.format(down),
n_estimators=len(FILE_NMBRS)*100,
)
pm[0].plot(xlim=(950,1000), ylim=(0,0.05), plot_poly_timestamp=410.)
Save the pulser model for the efficiency simulation.
with open(PATH_PULSER_MODEL, 'wb') as f:
pickle.dump(pm, f)
We plot the TPA equivalent spectrum to find the CPE factor.
x_ranges = [(0, 4), (0, 6)]
y_ranges = [(0, 200), (0, 500)]
for c, xr, yr in zip(H5_CHANNELS, x_ranges, y_ranges):
dh_stream.show_values(group='events', key='tpa_equivalent', bins=200, idx0=c, range=xr, yran=yr,
xlabel='TPA Equivalent (V)', ylabel='Counts', show=False)
# plt.axvline(x=peak_pos[c], color='black', linestyle='dotted', linewidth=2.5, zorder=70) # run the cell below first!
plt.show()
Take the peaks and calculate a CPE factor.
iron_peak = PEAK_ENERGY
tpa_equivalents = dh_stream.get('events', 'tpa_equivalent')
tpa_peak_cut = ai.cuts.LogicalCut(tpa_equivalents[0] > 1.6)
tpa_peak_cut.add_condition(tpa_equivalents[0] < 2.4)
tpa_peak_cut.add_condition(tpa_equivalents[1] < 0.9)
peak_pos = []
cpe_factor = []
for c in H5_CHANNELS:
peak_pos.append(np.median(tpa_equivalents[c, tpa_peak_cut.get_flag()]))
print('Peak position Channel {}: {} V_TPA'.format(c, peak_pos[-1]))
cpe_factor.append(iron_peak/peak_pos[-1])
print('CPE factor Channel {}: {} keV/V_TPA'.format(c, cpe_factor[-1]))
We can overwrite the calculated recoil energies with the new ones, with a precise CPE factor.
for c in H5_CHANNELS:
dh_stream.include_values(values=dh_stream.get('events', 'tpa_equivalent')[c]*cpe_factor[c], naming='recoil_energy', channel=c, delete_old=True)
dh_stream.include_values(values=dh_stream.get('events', 'tpa_equivalent_sigma')[c]*cpe_factor[c], naming='recoil_energy_sigma', channel=c, delete_old=True)
Quality Cuts¶
We apply some logical cuts. Before that we have a look at the main parameter distribution.
ph = dh_stream.get('events', 'pulse_height')
ydata = dh_stream.get('events', 'rise_time')
# ydata = dh_stream.get('events', 'decay_time')
xlims = [(0,1.5), (0,0.8)]
ylims = [(-1,5), (-1,10)]
# ylims = [(-1,40), (-1,25)]
for c, x, y in zip(H5_CHANNELS, xlims, ylims):
plt.close()
xlims, ylims, I = ai.styles.scatter_img(x_data=ph[c], y_data=ydata[c], ylims=y, xlims=x, height=100, width=200, alpha=0.5)
ax_extent = list(xlims)+list(ylims)
ai.styles.use_cait_style()
plt.imshow(I,
vmin=0,
vmax=1,
cmap=plt.get_cmap('viridis'),
interpolation='lanczos',
aspect='auto',
extent=ax_extent,
)
#plt.axhline(y=0.25, color='red', linewidth=2) # three event types
#plt.axhline(y=0.3, color='red', linewidth=2) # three event types
#plt.axvline(x=2.5, color='green', linewidth=2) # saturation
#plt.axvline(x=0.5, color='green', linewidth=2) # noise
plt.grid(alpha=0.2)
plt.title(' Channel {}'.format(c))
plt.xlabel('Pulse Height (V)')
plt.ylabel('Rise Time (ms)')
# plt.ylabel('Decay Time (ms)')
plt.show()
for c in H5_CHANNELS:
dh_stream.show_values(group='events', key='mainpar', bins=200, idx0=c, idx2=3, #range=(-5e-6, 5e-6),
xlabel='Sample Index', ylabel='Counts', title='Maximum Position Channel {}'.format(c))
ph_cut = ai.cuts.LogicalCut(initial_condition=np.abs(dh_stream.get('events', 'slope')[0,:]) < 0.2) # slope
ph_cut.add_condition(np.abs(dh_stream.get('events', 'slope')[1,:]) < 0.2) # slope
ph_cut.add_condition(dh_stream.get('events', 'pulse_height')[0,:,0] < 1.6) # ph
ph_cut.add_condition(dh_stream.get('events', 'pulse_height')[1,:,0] < 0.3) # ph
ph_cut.add_condition(dh_stream.get('events', 'onset')[0,:,3] < 50) # max pos phonon
ph_cut.add_condition(dh_stream.get('events', 'onset')[0,:,3] > 10) # max pos phonon
ph_cut.add_condition(dh_stream.get('events', 'pulse_height')[0,:] < MAXIMAL_EVENT_HEIGHTS[0]) # saturation
ph_cut.add_condition(dh_stream.get('events', 'of_ph')[0,:] > FIT_THRESHOLDS[0]) # threshold
print('Surviving Event Ratio: ', ph_cut.counts(), np.sum(ph_cut.get_flag())/len(ph_cut.get_flag()))
for c in H5_CHANNELS:
dh_stream.apply_logical_cut(cut_flag=ph_cut.get_flag(),
naming='ph_cut',
channel=c,
type='events',
delete_old=False)
Cut Efficiency¶
We determine the cut efficiency, that we need in the high level analysis to re-weight the event counts. This is the point at which you should run the script for the cut efficiency simulation.
dh_eff = ai.DataHandler(run=RUN,
module=MODULE,
channels=RDT_CHANNELS)
dh_eff.set_filepath(path_h5=PATH_PROC_DATA,
fname=FNAME_EFFICIENCY,
appendix=False)
dh_eff.calc_rate_cut(intervals=dh_stream.get('metainfo', 'rate_stable'))
nmbr_bins = 500
total_hours = 1000
good_intervals = dh_eff.get('metainfo', 'rate_stable')
dh_eff.show_values(group='events', key='hours', bins=nmbr_bins,
xlabel='Time (h)', ylabel=f'Counts/({total_hours/nmbr_bins*60:.1f} min)', title='Event Rate',
show=False)
for iv in good_intervals:
plt.axvspan(xmin=iv[0], xmax=iv[1], color='green', alpha=0.3, zorder=80)
# plt.xlim((400,450))
plt.show()
for c in H5_CHANNELS:
dh_eff.calc_controlpulse_stability(channel=c, instable_iv=dh_stream.get('metainfo', f'controlpulse_instable_ch{c}'))
surviving = ai.cuts.LogicalCut(initial_condition=np.abs(dh_eff.get('events', 'slope')[0]) < 0.1)
surviving.add_condition(np.abs(dh_eff.get('events', 'onset')[0]) < 20) # max pos phonon
surviving.add_condition(np.abs(dh_eff.get('events', 'rise_time')[0]) < 2.5) # rise time
surviving.add_condition(np.abs(dh_eff.get('events', 'rise_time')[0]) > 0.5) # rise time
surviving.add_condition(dh_eff.get('events', 'rate_cut')) # rate cut
surviving.add_condition(dh_eff.get('events', 'controlpulse_stability')[0]) # stability cut
surviving.add_condition(dh_eff.get('events', 'pulse_height')[0,:] < MAXIMAL_EVENT_HEIGHTS[0]) # saturation
surviving.add_condition(dh_eff.get('events', 'of_ph')[0,:] > FIT_THRESHOLDS[0]) # threshold
print('Surviving Event Ratio: ', surviving.counts(), np.sum(surviving.get_flag())/len(surviving.get_flag()))
for c in H5_CHANNELS:
dh_eff.apply_logical_cut(cut_flag=surviving.get_flag(),
naming='surviving',
channel=c,
type='events',
delete_old=False)
surviving = dh_eff.get('events', 'surviving')
with open(PATH_PULSER_MODEL, 'rb') as f:
pm = pickle.load(f)
# run this only if you havent during the efficiency loop
dh_eff.calc_calibration(starts_saturation=MAXIMAL_EVENT_HEIGHTS,
cpe_factor=CPE_FACTOR,
poly_order=3,
plot=False,
method='of',
pulser_models=pm,
name_appendix_energy='_reconstructed',
)
dh_eff.calc_calibration(starts_saturation=MAXIMAL_EVENT_HEIGHTS,
cpe_factor=CPE_FACTOR,
poly_order=3,
plot=False,
method='true_ph',
pulser_models=pm,
name_appendix_energy='_true',
)
for c in [H5_CHANNELS[0]]:
efficiency, counts, bins = dh_eff.show_efficiency(channel=c,
cut_flag=surviving[c],
which_quantity='recoil_energy_true',
bins=200,
title='Cut Efficiency',
xlabel='True Energy (keV)',
ylabel='Survival Probability',
range=[1e-3, 2.6],
show=True,
xscale='linear')
… and only for the low energy region:
for c in [H5_CHANNELS[0]]:
efficiency, counts, bins = dh_eff.show_efficiency(channel=c,
cut_flag=surviving[c],
which_quantity='recoil_energy_true',
bins=200,
title='Cut Efficiency',
xlabel='True Energy (keV)',
ylabel='Survival Probability',
show=True,
xscale='linear',
range=(1e-3, 0.5),)
Fit trigger efficiency.
for c in H5_CHANNELS:
ai.fit.fit_trigger_efficiency(binned_energies=bins, survived_fraction=efficiency,
a1_0=0.09, a0_0=1, a2_0=0.01, plot=True,
title='Trigger Efficiency Channel {}'.format(c),
xlim=(0,0.2),
)
Write the time dependent cut efficiency to a text file for limit calculations.
np.savetxt(PATH_PROC_DATA + 'xy_files/' + MODULE + '_cuteff.xy',
np.column_stack([dh_eff.get('events', 'hours')*60*60*1e6,
dh_eff.get('events', 'recoil_energy_true')[0],
dh_eff.get('events', 'recoil_energy_reconstructed')[0],
dh_eff.get('events', 'surviving')[0].astype(int),
(dh_eff.get('events', 'of_ph')[0] > FIT_THRESHOLDS[0]).astype(int)]))
Spectrum and Light Yield Plot¶
We plot an energy spectrum.
x_ranges = [(0, 8), (0, 80)]
y_ranges = [(0, 250), (0, 150)]
bins=200
for c, xr, yr in zip(H5_CHANNELS, x_ranges, y_ranges):
dh_stream.show_values(group='events', key='recoil_energy', bins=bins, idx0=c, cut_flag=ph_cut.get_flag(), range=xr, yran=yr,
xlabel='Recoil Energy (keV)', ylabel='Counts / {} eV'.format(int((xr[1]-xr[0])/bins*1000)),
save_path='plots/{}/{}_Spectrum_LE_Channel{}.pdf'.format(MODULE, MODULE, c))
… and only the low-energy part.
x_ranges = [(0, 0.5), (0, 80)]
y_ranges = [(0, 50), (0, 150)]
bins=150
for c, xr, yr in zip([0], x_ranges, y_ranges):
dh_stream.show_values(group='events', key='recoil_energy', bins=bins, idx0=c, cut_flag=ph_cut.get_flag(), range=xr, yran=yr,
xlabel='Recoil Energy (keV)', ylabel='Counts / {} eV'.format(int((xr[1]-xr[0])/bins*1000)),
save_path='plots/{}/{}_Spectrum_LE_Channel{}.pdf'.format(MODULE, MODULE, c))
Sometimes we see the peaks better, when we perform density estimation with a Gaussian kernal and the individual recoil energy uncertainties.
recoil_energy = dh_stream.get('events', 'recoil_energy')
recoil_energy = recoil_energy[:, ph_cut.get_flag()]
recoil_energy_sigma = dh_stream.get('events', 'recoil_energy_sigma')
recoil_energy_sigma = recoil_energy_sigma[:, ph_cut.get_flag()]
x_ranges = [(0, 10), (0, 100)]
y_ranges = [(0, 35000), (0, 3500)]
for c, xr, yr in zip(H5_CHANNELS, x_ranges, y_ranges):
# density estimation
x_d = np.linspace(xr[0], xr[1], 1000)
density = sum(norm(loc=m, scale=s).pdf(x_d) for m,s in zip(recoil_energy[c], recoil_energy_sigma[c]))/EXPOSURE
plt.close()
ai.styles.use_cait_style(dpi=150)
plt.fill_between(x_d, density, alpha=0.8, zorder=25, label='Density Estimate')
plt.plot(recoil_energy[c], np.full_like(recoil_energy[c], -0.1), '|k', markeredgewidth=1, label='Counts')
ai.styles.make_grid()
plt.xlim(xr)
plt.ylim(yr)
plt.xlabel('Recoil Energy (keV)')
plt.ylabel('Density (1/(keV kg days))')
plt.legend()
plt.savefig('plots/{}/{}_Density_Channel{}.pdf'.format(MODULE, MODULE, c))
plt.show()
We single out the iron lines.
recoil_energy = dh_stream.get('events', 'recoil_energy')
recoil_energy_sigma = dh_stream.get('events', 'recoil_energy_sigma')
iron_lower = ai.cuts.LogicalCut(recoil_energy[0] > 5)
iron_lower.add_condition(recoil_energy[0] < 6.2)
iron_lower.add_condition(ph_cut.get_flag())
iron_upper = ai.cuts.LogicalCut(recoil_energy[0] > 6.2)
iron_upper.add_condition(recoil_energy[0] < 6.7)
iron_upper.add_condition(ph_cut.get_flag())
rec_lower = recoil_energy[:, iron_lower.get_flag()]
rec_sigma_lower = recoil_energy_sigma[:, iron_lower.get_flag()]
rec_upper = recoil_energy[:, iron_upper.get_flag()]
rec_sigma_upper = recoil_energy_sigma[:, iron_upper.get_flag()]
x_ranges = [(5, 7), (0, 20)]
y_ranges = [(0, 30000), (0, 2000)]
for c, xr, yr in zip(H5_CHANNELS, x_ranges, y_ranges):
print('Mean lower line, channel {}: {}'.format(c, np.mean(rec_lower[c])))
print('Mean upper line, channel {}: {}'.format(c, np.mean(rec_upper[c])))
# density estimation lower iron
x_d = np.linspace(xr[0], xr[1], 1000)
density_lower = sum(norm(loc=m, scale=s).pdf(x_d) for m,s in zip(rec_lower[c], rec_sigma_lower[c]))/EXPOSURE
# density estimation upper iron
density_upper = sum(norm(loc=m, scale=s).pdf(x_d) for m,s in zip(rec_upper[c], rec_sigma_upper[c]))/EXPOSURE
plt.close()
ai.styles.use_cait_style(dpi=150)
# lower iron
plt.fill_between(x_d, density_lower, alpha=0.8, zorder=25, label='K alpha')
plt.plot(rec_lower[c], np.full_like(rec_lower[c], -0.1), '|k', markeredgewidth=1)
# upper iron
plt.fill_between(x_d, density_upper, alpha=0.8, zorder=25, label='K beta')
plt.plot(rec_upper[c], np.full_like(rec_upper[c], -0.1), '|k', markeredgewidth=1)
ai.styles.make_grid()
plt.xlim(xr)
plt.ylim(yr)
plt.xlabel('Recoil Energy (keV)')
plt.ylabel('Density (1/(keV kg days))')
plt.legend()
plt.show()
And we do a light yield plot.
recoil_energy_corr = dh_stream.get('events', 'recoil_energy')
recoil_energy_corr = recoil_energy_corr[:, ph_cut.get_flag()]
plt.close()
xlims, ylims, I = ai.styles.scatter_img(x_data=recoil_energy_corr[0],
y_data=recoil_energy_corr[1]/recoil_energy_corr[0],
ylims=(-4, 10),
xlims=(0, 10),
height=200,
width=250,
alpha=0.9)
ax_extent = list(xlims)+list(ylims)
ai.styles.use_cait_style()
plt.imshow(I,
vmin=0,
vmax=1,
cmap=plt.get_cmap('viridis'),
interpolation='lanczos',
aspect='auto',
extent=ax_extent,
)
plt.grid(alpha=0.2)
plt.title('Light Yield Plot')
plt.xlabel('Recoil Energy (keV)')
plt.ylabel('Light Yield (a.u.)')
#plt.colorbar()
plt.savefig('plots/{}/{}_LY.pdf'.format(MODULE, MODULE))
plt.show()
What is the event rate of low energetic events?
only_le = ai.cuts.LogicalCut(initial_condition=dh_stream.get('events', 'recoil_energy')[0] < 0.5) # cut at 0.5 keV
only_le.add_condition(ph_cut.get_flag())
print('Total Counts: ', ph_cut.counts())
print('Counts below 0.5 keV: ', only_le.counts())
dh_stream.show_values(group='events', key='hours', bins=100,
xlabel='Time (h)', ylabel='Counts', title='Event Rate below 0.5 keV', cut_flag=only_le.get_flag())
Write the spectrum to a xy file to process it for limit calculations.
recoil_energy_corr = dh_stream.get('events', 'recoil_energy')
recoil_energy_corr = recoil_energy_corr[:, ph_cut.get_flag()]
data = [recoil_energy_corr[0],
recoil_energy_corr[1]/recoil_energy_corr[0]]
ai.data.write_xy_file(filepath=PATH_PROC_DATA + 'xy_files/' + MODULE + '_spectrum.xy',
data=data,
title='Run ' + RUN + ' ' + MODULE + ' Spectrum',
axis=['Recoil Energy (keV)',
'Light Yield (a.u.)'])
data.append(dh_stream.get('events', 'hours')[ph_cut.get_flag()])
ai.data.write_xy_file(filepath=PATH_PROC_DATA + 'xy_files/' + MODULE + '_spectrum_timedependent.xy',
data=data,
title='Run ' + RUN + ' ' + MODULE + 'Spectrum Time Dependent',
axis=['Recoil Energy (keV)',
'Light Yield (a.u.)',
'Time (h)'])
Write the times at which the individual files start and stop to an xy file.
dh_stream.generate_startstop()
ai.data.write_xy_file(filepath=PATH_PROC_DATA + 'xy_files/' + MODULE + '_startstop.xy',
data=dh_stream.get('metainfo', 'startstop_hours').T,
title='Run ' + RUN + ' ' + MODULE + ' Startstop',
axis=['Start Files (hours)',
'Stop Files (hours)'])
np.savetxt(PATH_PROC_DATA + 'xy_files/' + MODULE + '_startstop_tstamp.xy',
np.column_stack([dh_stream.get('metainfo', 'startstop_s')[:,0].astype(int),
dh_stream.get('metainfo', 'startstop_mus')[:,0].astype(int),
dh_stream.get('metainfo', 'startstop_s')[:,1].astype(int),
dh_stream.get('metainfo', 'startstop_mus')[:,1].astype(int)]))
Finally, what is the exposure?
dh_stream.exposure(detector_mass=DETECTOR_MASS, exclude_instable=True)
Finished.