Select Git revision
evaluation.cpython-311.pyc
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
PeakOTron.py 94.86 KiB
from HelperFunctions import GP_gain, GP_lbda, GP_mu, HistMean, HistStd, HistSkew
from HelperFunctions import LatexFormat, Linear
from AdditionalPDFs import gpoisson
from LossFunctions import *
from Model import DRM, Beta, Beta_Inf
from iminuit import Minuit
from scipy.interpolate import interp1d, UnivariateSpline, InterpolatedUnivariateSpline
from scipy.stats import norm
from astropy.stats import knuth_bin_width, freedman_bin_width, scott_bin_width
from matplotlib import cm
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.ticker import ScalarFormatter
from itertools import chain
from copy import deepcopy
import time
import matplotlib.pyplot as plt
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
########################################################################################################
## PEAKOTRON
########################################################################################################
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
class PeakOTron:
########################################################################################################
## __init__
########################################################################################################
##
## VARIABLES:
## n_call_minuit: number of calls to Minuit for each trial
## n_iterations_minuit: number of trial iterations if convergence is not achieved
## n_bootstrap: number of bootstraps for each statistic
## plot_figsize: figure size in inches for prefit figures
## plot_fontsize: fontsize for output plots
## plot_cmap: colour map for output plots
##
## This function defines the basic information stored by PeakOTron and the default fit value dictionary that is updated with the input parameters of ## the user.
##
########################################################################################################
def __init__(self,
verbose=False,
n_call_minuit=10000,
n_iterations_minuit = 10,
plot_figsize=(10,10),
plot_fontsize=25,
plot_legendfontsize=18,
plot_cmap = 'turbo',
eps = 1e-6
):
self._default_hist_data={
"count":None,
"density":None,
"bin_numbers":None,
"bin_centres":None,
"counts_bg":None,
"counts_bgsub":None,
"truncate_nsigma0_do_min":None,
"bw":None,
"peaks":{
"peak_position":None,
"peak_height":None,
"peak_index":None,
"peak_mean":None,
"peak_mean_error":None,
"peak_variance":None,
"peak_variance_error":None,
"peak_std_deviation":None,
"peak_std_deviation_error":None
}
}
self._default_fit_dict={
"Q_0":None,
"G":None,
"lbda":None,
"mu":None,
"sigma_0":None,
"sigma_1":None,
"DCR":None,
"tau":None,
"t_gate":None,
"t_0":None,
"k_max":None,
"len_pad":None
}
self._eps = eps
self._plot_figsize= plot_figsize
self._plot_fontsize= plot_fontsize
self._plot_legendfontsize= plot_legendfontsize
self._cmap = cm.get_cmap(plot_cmap)
self._len_pad = int(100)
self._verbose=verbose
self._n_call_minuit = n_call_minuit
self._n_iterations_minuit = n_iterations_minuit
self._max_n_peaks = None
self._max_n_peaks_dcr = 6
self._min_n_events_peak=100
self._m_DRM = None
self._time_prefit = None
self._time_fit = None
self._default_fit_kwargs={
"tau":None,
"tau_err":None,
"t_0":None,
"t_0_err":None,
"tau_R":None,
"tau_R_err":None,
"t_gate":None,
"t_gate_err":None,
"bw":None,
"bin_0":None,
"truncate_nsigma0_up":None,
"truncate_nsigma0_do":None,
"truncate_nsigma0_chi2scanlim":None,
"min_n_events_peak":None,
"prefit_only":False,
"bin_method":"knuth"
}
self.InitKwargs()
########################################################################################################
## InitKwargs()
########################################################################################################
##
## This function creates the fit dictionary and initialises the output dictionaries of the class.
##
########################################################################################################
def InitKwargs(self):
self._fit_kwargs = self._default_fit_kwargs.copy()
self._prefit_values = self._default_fit_dict.copy()
self._prefit_errors = self._default_fit_dict.copy()
self._fit_values= self._default_fit_dict.copy()
self._fit_errors = self._default_fit_dict.copy()
self._prefit_values_bin = self._default_fit_dict.copy()
self._prefit_errors_bin = self._default_fit_dict.copy()
self._fit_values_bin= self._default_fit_dict.copy()
self._fit_errors_bin = self._default_fit_dict.copy()
self._hist_data = self._default_hist_data.copy()
########################################################################################################
## SetMaxNPeaks(max_n_peaks)
########################################################################################################
##
## Typically, it is sufficient to determine the appropriate number of peaks to fit from mu.
## However, to make sure we don't accidentally create a fit with a large number of peaks,
## the default threshold of 20 peaks may be set with this function
##
########################################################################################################
def SetMaxNPeaks(self, max_n_peaks):
max_n_peaks = int(max_n_peaks)
if(max_n_peaks<4):
raise Exception("Maximum number of light peaks must be greater than 3.")
self._max_n_peaks = max_n_peaks
if(self._verbose):
print("Set maximum number of light peaks to {:d}.".format(self._max_n_peaks))
########################################################################################################
## SetMaxNPeaksDCR(max_n_peaks)
########################################################################################################
##
## Typically, it is sufficient to determine the appropriate number of peaks to fit from mu.
## However, to make sure we don't accidentally create a fit with a large number of peaks,
## the default threshold of 20 peaks may be set with this function
##
########################################################################################################
def SetMaxNPeaksDCR(self, max_n_peaks_dcr):
max_n_peaks_dcr = int(max_n_peaks_dcr)
if(max_n_peaks_dcr<4):
raise Exception("Maximum number of dark peaks must be greater than 3.")
self._max_n_peaks_dcr = max_n_peaks_dcr
if(self._verbose):
print("Set maximum number of dark peaks to {:d}.".format(self._max_n_peaks_dcr))
########################################################################################################
## GetMinuit()
########################################################################################################
##
## Returns Minuit object after fit. Note this object has fit results stored in bins, not in the unit parameter of the input.
##
########################################################################################################
def GetMinuit(self):
if(self._m_DRM is None):
raise Exception ("Fit has not yet been performed. Please first perform a fit.")
else:
return self._m_DRM
########################################################################################################
## GetFitTime()
########################################################################################################
##
## Returns fit time, for prefit and fit.
##
########################################################################################################
def GetFitTime(self, prefit=False):
if(prefit):
if(self._time_prefit is None):
raise Exception ("Fit has not yet been performed. Please first perform a fit.")
else:
return self._time_prefit
else:
if(self._time_fit is None):
raise Exception ("Fit has not yet been performed. Please first perform a fit.")
else:
return self._time_fit
########################################################################################################
## IsValid()
########################################################################################################
##
## Returns whether or not Minuit returned valid.
##
########################################################################################################
def IsValid(self):
if(self._m_DRM is None):
raise Exception ("Fit has not yet been performed. Please first perform a fit.")
else:
return self._m_DRM.valid
########################################################################################################
## GetFitResults()
########################################################################################################
##
## Returns fit results and errors in the input charge unit.
##
########################################################################################################
def GetFitResults(self, bin_units=True):
if(self._m_DRM is None):
raise Exception ("Fit has not yet been performed. Please first perform a fit.")
else:
if(bin_units):
temp_values = deepcopy(self._fit_values_bin)
temp_errors = deepcopy(self._fit_errors_bin)
else:
temp_values = deepcopy(self._fit_values)
temp_errors = deepcopy(self._fit_errors)
## This part rescales the afterpulse to an infinite gate length.
## Please note that the effectiveness of the scaling factor depends on the ability to determine tau_Ap.
## This may be challenging for some spectra where there is insufficient information between the peaks to calculate tau_Ap.
## In that case, it is suggested that an appropriate tau_Ap is assumed (i.e. from overvoltages allowing an accurate assessment of tau_Ap)
## so that this factor can be reconstructed.
p_Ap_inf_sf = Beta_Inf(temp_values["tau_R"], temp_values["tau_Ap"])/Beta(temp_values["tau_R"], temp_values["tau_Ap"], temp_values["t_gate"])
temp_values["p_Ap_inf"] = temp_values["p_Ap"]*p_Ap_inf_sf
temp_errors["p_Ap_inf"] = temp_errors["p_Ap"]*p_Ap_inf_sf ## TO DO: propogate errors correctly.
t_fit = self.GetFitTime(prefit=False)
t_prefit = self.GetFitTime(prefit=True)
temp_values["fit_time"] = t_fit
temp_values["prefit_time"] = t_prefit
return temp_values, temp_errors
########################################################################################################
## GetPrefitResults()
########################################################################################################
##
## Returns pre-fit results and errors in the input charge unit.
##
########################################################################################################
def GetPrefitResults(self, bin_units=True):
if(bin_units):
return self._prefit_values_bin, self._prefit_errors_bin
else:
return self._prefit_values, self._prefit_errors
########################################################################################################
## GetPeakInfo()
########################################################################################################
##
## Returns peak information obtained during the fitting procedure.
##
########################################################################################################
def GetPeakInfo(self):
if(self._m_DRM is None):
raise Exception ("Fit has not yet been performed. Please first perform a fit.")
else:
return self._hist_data
########################################################################################################
## GetChi2(prefit=False)
########################################################################################################
##
## VARIABLES:
## prefit(False): Evaluate the model using the prefit or fitted values
##
##
##
## Calculates the Chi2 of the resultant fit relative to the original density estimate.
########################################################################################################
def GetChi2(self, pars=None, prefit=False, n_sigma_lims=[None, None]):
_x_all = self._hist_data["bin_numbers"]
_y_all = self._hist_data["counts"]
_y_err_all = self._hist_data["counts_error"]
_Q_0_prefit = self._prefit_values_bin["Q_0"]
_sigma_0_prefit = self._prefit_values_bin["sigma_0"]
if(pars is None):
_y_hat_all = self.GetModel(_x_all, prefit=prefit, bin_units=True)
else:
_y_hat_all = DRM(_x_all, **pars)
cond_count = (_y_all>1)
if(n_sigma_lims[0] is None):
if(self._hist_data["truncate_nsigma0_do_min"] is not None):
_lim_nsigma0_low = _Q_0_prefit -self._hist_data["truncate_nsigma0_do_min"]*_sigma_0_prefit
else:
_lim_nsigma0_low = -np.inf
else:
_lim_nsigma0_low = _Q_0_prefit - n_sigma_lims[0]*_sigma_0_prefit
if(n_sigma_lims[1] is None):
_lim_nsigma0_hi = np.inf
else:
_lim_nsigma0_hi = _Q_0_prefit + n_sigma_lims[1]*_sigma_0_prefit
cond_nsigma0_low = (_x_all > _lim_nsigma0_low)
cond_nsigma0_hi = (_x_all < _lim_nsigma0_hi)
conds = cond_count & cond_nsigma0_low & cond_nsigma0_hi
_x = _x_all[conds]
_y = _y_all[conds]
_y_err = _y_err_all[conds]
_y_hat = _y_hat_all[conds]
_y_hat_err = np.sqrt(_y_hat_all[conds])
chi2 = np.sum(((_y_hat - _y)/_y_hat_err)**2)
ndof = len(_x)-10
return chi2, ndof
########################################################################################################
## GetModel(x, prefit=False)
########################################################################################################
##
## VARIABLES:
## x: the x values over which to evaluate the model.
## prefit (False): whether or not to use the prefit values in the model evaluation.
##
##
## Evaluates model according to the fitted parameters and returns the Detector Response Model
########################################################################################################
def GetModel(self, x, prefit=False, bin_units=False):
if(self._m_DRM is None):
raise Exception ("Fit has not yet been performed. Please first perform a fit.")
else:
if(prefit):
if(bin_units):
return DRM(x, **self._prefit_values_bin)
else:
return DRM(x, **self._prefit_values)
else:
if(bin_units):
return DRM(x, **self._fit_values_bin)
else:
return DRM(x, **self._fit_values)
########################################################################################################
## GetBins(data, bw, bin_method)
########################################################################################################
##
## VARIABLES:
## data: the input data, either a 1D array of individual discharges or a 2D array of [bin centres, counts]
## bw: the bin width. If this is None, then the bin_method will be used to choose the binning if 1D
## bin_method: the method for binning 1D data. This can be 'knuth', 'freedman', or 'scott'
##
## Automated calculation of bins and bin widths for 2D prebinned charge measurements or 1D charge measurements.
## Note: because spectra from scopes have been observed with variable bin widths, the program also rebins 2D histograms using the same bin width as the original spectrum.
########################################################################################################
def GetBins(self, data, bw, bin_method, bin_0):
if(data.ndim == 1):
if(self._verbose):
print("Data is assumed 1D. Assuming list of charges.")
if(bw is None):
if(self._verbose):
print("Bin width not set. Binning with {:s}".format(bin_method))
if(bin_method == "knuth"):
bw = knuth_bin_width(data)/20
elif(bin_method == "freedman"):
bw = freedman_bin_width(data)
elif(bin_method == "scott"):
bw = scott_bin_width(data)
else:
raise Exception("Binning method not recognised. Please select bin_method from 'knuth', 'freedman' or 'scott'")
x_min, x_max = data.min(), data.max()
elif(data.ndim == 2):
if(self._verbose):
print("Data is assumed 2D. Assuming input is a histogram, with columns [bin_centre, counts].")
_x, _y = data[:,0], data[:,1]
data = np.array(list(chain.from_iterable([[__x]*int(__y) for __x,__y in zip(_x,_y)])))
f_bin = interp1d(np.arange(0, len(_x)), _x)
if(bw is None):
bw = abs(f_bin(1) - f_bin(0))
idx_nonzero = np.squeeze(np.argwhere((_y>0)))
x_min, x_max = _x[np.min(idx_nonzero)]-bw/2, _x[np.max(idx_nonzero)]+bw/2
else:
raise Exception("Data is neither 1D nor 2D. Please provide the correct dimensions of data to the fit program.")
if(bin_0 is not None):
if(self._verbose):
print("Setting first bin to be {:3.3f} input charge units".format(bin_0))
x_min = bin_0
nbins = np.ceil((x_max - x_min) / bw)
nbins = max(1, nbins)
bins = x_min + bw * np.arange(nbins + 1)
return data, bw, nbins, bins
def PlotHistogram(self, figsize=(10.0,10.0), save_directory=None, label="Histogram", x_label="ADC"):
fig = plt.figure(figsize=figsize)
ax = plt.subplot(111)
counts = self._hist_data["counts"]
bin_numbers = self._hist_data["bin_numbers"]
bin_centres = self._hist_data["bin_centres"]
ax.step(bin_numbers,
counts,
#label="{:s} \n $N_{{events}}$ = {:s}".format(label, LatexFormat(float(np.sum(counts)))),
label = "{:s}".format(label),
color="C0",
where="mid",
lw=5)
ax.grid(which="both")
ax.set_ylabel("Counts", fontsize=self._plot_fontsize)
ax.set_xlabel("Bin", fontsize=self._plot_fontsize)
ax.tick_params(axis="x", labelsize=self._plot_fontsize, rotation=45)
ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.set_ylim([1,None])
ax.set_yscale("log")
ax.legend(fontsize=self._plot_legendfontsize)
_f = interp1d(bin_numbers, bin_centres, bounds_error=False, fill_value="extrapolate")
_f_inv = interp1d(bin_centres, bin_centres, bounds_error=False, fill_value="extrapolate")
# secax = ax.secondary_xaxis('top', functions=(_f, _f_inv))
# secax.set_xlabel('{:s}'.format(x_label), fontsize=self._plot_fontsize)
# secax.tick_params(axis="x", labelsize=self._plot_fontsize, rotation=45)
# mf = ScalarFormatter(useMathText=True)
# mf.set_powerlimits((-2,2))
# secax.xaxis.set_major_formatter(mf)
# offset = secax.xaxis.get_offset_text()
# offset.set_size(int(0.8*self._plot_fontsize))
if(save_directory is not None):
print("Saving figure to {0}...".format(save_directory))
fig.savefig(save_directory)
if(display):
plt.pause(0.01)
fig.show()
else:
plt.close(fig)
def PlotFFT(self, figsize=(10.0,10.0), save_directory=None, label="Histogram"):
fig = plt.figure(figsize=figsize)
ax = plt.subplot(111)
fft_freq = self._hist_data["fft_freq"]
counts = self._hist_data["counts"]
fft_amplitude = self._hist_data["fft_amplitude"]
G_fft = self._hist_data["G_fft"]
inv_G_fft = 1/G_fft
ax.step(fft_freq,
fft_amplitude,
color="C0",
where="mid",
#label="{:s} \n $N_{{events}}$ = {:s}".format(label, LatexFormat(float(np.sum(counts)))),
label = "{:s}".format(label),
lw=5)
ax.axvline(x=inv_G_fft, color="purple", linestyle="--", lw=4, label="$G^{{*}}_{{FFT}}$ = {:3.3f} Bin".format(G_fft))
ax.set_yscale("log")
ax.set_xscale("log")
ax.grid(which="both")
ax.set_xlabel("Bin in Fourier Space", fontsize=25)
ax.set_ylabel("Power Spectral Density", fontsize=25)
ax.tick_params(axis="x", labelsize=self._plot_fontsize, rotation=45)
ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.legend(fontsize=self._plot_legendfontsize)
if(save_directory is not None):
print("Saving figure to {0}...".format(save_directory))
fig.savefig(save_directory)
if(display):
plt.pause(0.01)
fig.show()
else:
plt.close(fig)
def PlotBGSub(self, figsize=(10,10), save_directory=None, plot_bgsub=False, label="Histogram"):
fig = plt.figure(figsize=figsize)
ax = plt.subplot(111)
counts = self._hist_data["counts"]
counts_bg = self._hist_data["counts_bg"]
counts_bgsub = self._hist_data["counts_bgsub"]
bin_numbers = self._hist_data["bin_numbers"]
if(plot_bgsub):
ax.step(bin_numbers,
counts,
#label="{:s} \n $N_{{events}}$ = {:s}".format(label, LatexFormat(float(np.sum(counts)))),
label = "{:s}".format(label),
color="C0",
where="mid",
lw=5)
ax.step(bin_numbers,
counts_bgsub,
#label="{:s} - BG \n $N_{{events}}$ = {:s}".format(label, LatexFormat(float(np.sum(counts_bgsub)))),
label = "{:s} - Background".format(label),
color="r",
where="mid",
lw=5)
else:
ax.step(bin_numbers,
counts,
#label="{:s} \n $N_{{events}}$ = {:s}".format(label,LatexFormat(float(np.sum(counts)))),
label = "{:s}".format(label),
color="C0",
where="mid",
lw=5)
ax.fill_between(bin_numbers[counts_bg>0],
counts_bg[counts_bg>0],
#label="{:s}, BG \n $N_{{events}}$ = {:s}".format(label, LatexFormat(float(np.sum(counts_bg)))),
label = "{:s}, Background".format(label),
color='none',
hatch="///",
edgecolor="r",
step="mid",
lw=5)
ax.grid(which="both")
ax.set_xlabel("Bin", fontsize=25)
ax.set_ylabel("Counts", fontsize=self._plot_fontsize)
ax.tick_params(axis="x", labelsize=self._plot_fontsize, rotation=45)
ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.set_ylim([1,None])
ax.set_yscale("log")
ax.legend(fontsize=self._plot_legendfontsize)
if(save_directory is not None):
print("Saving figure to {0}...".format(save_directory))
fig.savefig(save_directory)
if(display):
plt.pause(0.01)
fig.show()
else:
plt.close(fig)
def PlotPeaks(self, figsize=(10.0,10.0), save_directory=None, label="Histogram"):
fig = plt.figure(figsize=figsize)
ax = plt.subplot(111)
counts = self._hist_data["counts"]
counts_bgsub = self._hist_data["counts_bgsub"]
bin_numbers = self._hist_data["bin_numbers"]
Q_0_est = self._hist_data["Q_0_est"]
peak_position = self._hist_data["peaks"]["peak_position"]
ax.step(bin_numbers,
counts,
#label="{:s} \n $N_{{events}}$ = {:s}".format(label, LatexFormat(float(np.sum(counts)))),
label = "{:s}".format(label),
color="C0",
lw=5,
where="mid",
zorder=-5)
ax.step(bin_numbers, counts_bgsub,
#label="{:s} - BG \n $N_{{events}}$ = {:s}".format(label, LatexFormat(float(np.sum(counts_bgsub)))),
label="{:s} - Background".format(label),
color="r",
lw=5,
where="mid",
zorder=-4)
color_vals = [self._cmap(_v) for _v in np.linspace(0,0.75, len(peak_position))]
for peak_i, pp in enumerate(peak_position):
if(peak_i == 0):
annotation_text = "$Q_{{0}}$"
color = "purple"
ls="--"
ax.axvline(x=pp,
color=color,
linestyle=ls,
lw=2.5,
label = '{:s}'.format(annotation_text)
)
else:
annotation_text = "Peak Position {:d}".format(peak_i)
color = color_vals[peak_i]
ls="-"
ax.axvline(x=pp,
color=color,
linestyle=ls,
lw=2.5
)
# plt.text(10.1,0, annotation_text, color=color, fontsize=self._plot_legend,rotation=90)
ax.scatter([Q_0_est],
[counts[np.argmin(abs(bin_numbers-Q_0_est))]],
marker="v",
s=200,
color="purple",
zorder=5,
label= "$Q^{est}_{{0}}$")
ax.set_xlabel("Bin", fontsize=self._plot_fontsize)
ax.set_ylabel("Counts", fontsize=self._plot_fontsize)
ax.tick_params(axis="x", labelsize=self._plot_fontsize, rotation=45)
ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.grid(which="both")
ax.legend(fontsize=self._plot_legendfontsize)
if(save_directory is not None):
print("Saving figure to {0}...".format(save_directory))
fig.savefig(save_directory)
if(display):
plt.pause(0.01)
fig.show()
else:
plt.close(fig)
def PlotVarianceFit(self, figsize=(10.0, 10.0), save_directory=None, y_limits=[None, None], x_limits=[None, None]):
fig = plt.figure(figsize=figsize)
ax = plt.subplot(111)
peak_index = self._hist_data["peaks"]["peak_index"]
peak_variance = self._hist_data["peaks"]["peak_variance"]
peak_variance_error = self._hist_data["peaks"]["peak_variance_error"]
sigma_0 = self._prefit_values_bin["sigma_0"]
dsigma_0 = self._prefit_errors_bin["sigma_0"]
sigma_1 = self._prefit_values_bin["sigma_1"]
dsigma_1 = self._prefit_errors_bin["sigma_1"]
ax.errorbar( peak_index,
peak_variance,
yerr=peak_variance_error,
fmt="o",
lw=3,
markersize=10,
color="C0",
label="Variances of Peaks"
)
ax.plot(peak_index,
Linear(peak_index,
sigma_1**2,
sigma_0**2
),
color="C0",
linestyle="--",
lw=3,
label = "Linear Fit: \n $\sigma^{{2}}_{{0}}$ = {:s} $\pm$ {:s} Bin$^{{2}}$ \n $\sigma^{{2}}_{{1}}$ = {:s} $\pm$ {:s} Bin$^{{2}}$".format(
LatexFormat(sigma_0**2),
LatexFormat(2*sigma_0*dsigma_0),
LatexFormat(sigma_1**1),
LatexFormat(2*sigma_0*dsigma_1),
),
)
ax.set_xlabel("Peak Index", fontsize=self._plot_fontsize)
ax.set_ylabel("$\sigma^{2}_{peak}$ [Bin$^{2}$]", fontsize=self._plot_fontsize)
ax.tick_params(axis="x", labelsize=self._plot_fontsize, rotation=45)
ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.grid(which="both")
ax.set_ylim(y_limits)
ax.set_xlim(x_limits)
ax.legend(fontsize=self._plot_legendfontsize)
if(save_directory is not None):
print("Saving figure to {0}...".format(save_directory))
fig.savefig(save_directory)
if(display):
plt.pause(0.01)
fig.show()
else:
plt.close(fig)
def PlotMeanFit(self, figsize=(10.0, 10.0), y_limits=[None, None], x_limits=[None, None], save_directory=None):
fig = plt.figure(figsize=figsize)
ax = plt.subplot(111)
bin_numbers = self._hist_data["bin_numbers"]
peak_index = self._hist_data["peaks"]["peak_index"]
peak_mean = self._hist_data["peaks"]["peak_mean"]
peak_mean_error = self._hist_data["peaks"]["peak_mean_error"]
Q_0 = self._prefit_values_bin["Q_0"]
dQ_0 = self._prefit_errors_bin["Q_0"]
G = self._hist_data["G_fft"]
ax.errorbar( peak_index,
peak_mean,
yerr=peak_mean_error,
fmt="o",
markersize=10,
lw=3,
color="C0",
label="Means of Peaks"
)
ax.plot( peak_index,
Linear(peak_index,
G,
Q_0,
),
color="C0",
linestyle="--",
lw=3,
markersize=10,
label = "Linear Fit: \n $Q_{{0}}$ = {:s} $\pm$ {:s} Bin \n $G^{{*}}_{{FFT}}$ = {:s} Bin".format(
LatexFormat(Q_0),
LatexFormat(dQ_0),
LatexFormat(G),
),
)
ax.set_xlabel("Peak Index", fontsize=self._plot_fontsize)
ax.set_ylabel("Peak Position [Bin]", fontsize=self._plot_fontsize)
# ax.set_ylabel("$\\langle Q_{{peak}} \\rangle$", fontsize=self._plot_fontsize)
ax.tick_params(axis="x", labelsize=self._plot_fontsize, rotation=45)
ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.grid(which="both")
ax.set_ylim(y_limits)
ax.set_xlim(x_limits)
ax.legend(fontsize=self._plot_legendfontsize)
if(save_directory is not None):
print("Saving figure to {0}...".format(save_directory))
fig.savefig(save_directory)
if(display):
plt.pause(0.01)
fig.show()
else:
plt.close(fig)
def PlotDCR(self, figsize=(10.0, 10.0), save_directory=None, label = "Histogram"):
fig = plt.figure(figsize=figsize)
ax = plt.subplot(111)
counts = self._hist_data["counts"]
bin_numbers = self._hist_data["bin_numbers"]
bw = self._hist_data["bw"]
Q_0 = self._prefit_values_bin["Q_0"]
G = self._prefit_values_bin["G"]
DCR = self._prefit_values_bin["DCR"]
dDCR = self._prefit_errors_bin["DCR"]
K = (bin_numbers - Q_0)/G
cond_NN = (K>=0.45)&(K<=0.55)
cond_Nc = (K<=0.5)
cond_DCR = (K<=2.5)
NN_sum = max(np.sum(counts[cond_NN]),1)
NN = np.mean(counts[cond_NN])
if(NN!=NN):
NN = float(counts[np.argmin(abs(K-0.5))])
NN_sum = float(NN)
Nc = max(float(np.sum(counts[cond_Nc])), 1)
ax.step(K[cond_DCR],
counts[cond_DCR],
label="{:s}".format(label),
color="C0",
where="mid",
lw=5)
ax.fill_between(
K[cond_Nc],
counts[cond_Nc],
color='none',
hatch="///",
edgecolor="r",
lw=3,
step="mid",
#label = "$N_{{K \leq 0.5}}$={:d} events".format(int(Nc)) )
label = "$N_{{K \leq 0.5}}$")
ax.fill_between(
K[cond_NN],
counts[cond_NN],
color='green',
edgecolor="g",
alpha=0.3,
step="mid",
#label = "$\\langle N \\rangle_{{0.45 < K \leq 0.55}}$={:3.1f} events".format(NN))
label = "$\\langle N \\rangle_{{0.45 < K \leq 0.55}}$")
legend_title="$DCR$={:s} $\pm$ {:s} MHz".format(LatexFormat(DCR*1e3), LatexFormat(dDCR*1e3))
ax.set_xlabel("K [p.e.]", fontsize=self._plot_fontsize)
ax.grid(which="both")
ax.set_ylabel("Counts", fontsize=self._plot_fontsize)
ax.tick_params(axis="x", which='both', labelsize=self._plot_fontsize, rotation=45)
ax.tick_params(axis="y", which='both', labelsize=self._plot_fontsize)
ax.set_yscale("log")
ax.legend(title=legend_title, title_fontsize=int(self._plot_legendfontsize), fontsize=int(self._plot_legendfontsize))
if(save_directory is not None):
print("Saving figure to {0}...".format(save_directory))
fig.savefig(save_directory)
if(display):
plt.pause(0.01)
fig.show()
else:
plt.close(fig)
def PlotFit(self,
figsize=(10.0,10.0),
labelsize=None,
ticksize=None,
titlesize=None,
legendsize=None,
title=None,
xlabel = None,
scaled=False,
save_directory=None,
label=None,
y_limits = [1, None],
residual_limits = [-5, 5],
residualticksize=None,
linewidth_main = 5,
x_limits = [None, None],
display=True,
prefit=False,
plot_in_bins=False
):
if(labelsize is None):
labelsize=self._plot_fontsize
if(ticksize is None):
ticksize=self._plot_fontsize
if(titlesize is None):
titlesize = self._plot_fontsize
if(legendsize is None):
legendsize = self._plot_legendfontsize
if(residualticksize is None):
residualticksize = int(0.8*self._plot_fontsize)
if(plot_in_bins):
x_all = self._hist_data["bin_numbers"]
bw = 1
else:
x_all = self._hist_data["bin_centres"]
bw = self._hist_data["bw"]
y_all = self._hist_data["counts"]
y_err_all = self._hist_data["counts_error"]
Q_0_prefit = self._prefit_values_bin["Q_0"]
sigma_0_prefit = self._prefit_values_bin["sigma_0"]
if(self._hist_data["truncate_nsigma0_do_min"] is not None):
lim_nsigma0 = Q_0_prefit - self._hist_data["truncate_nsigma0_do_min"]*sigma_0_prefit
cond_sel_nsigma0 = (self._hist_data["bin_numbers"]>lim_nsigma0)
else:
cond_sel_nsigma0 = np.full(len(x_all), True).astype(bool)
x = x_all[cond_sel_nsigma0]
y = y_all[cond_sel_nsigma0]
y_err = y_err_all[cond_sel_nsigma0]
y_hat = self.GetModel(x, prefit=prefit, bin_units=plot_in_bins)*bw
y_hat_err = np.sqrt(y_hat)
chi2, ndf = self.GetChi2(prefit=prefit)
if(xlabel is None):
if(plot_in_bins):
xlabel="Bin"
else:
xlabel = "Q"
fig = plt.figure(figsize=figsize)
gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])
ax0 = plt.subplot(gs[0])
if(label is None):
label = "Histogram"
# label="{:s},\n $N_{{events}}$ = {:s}".format(label, LatexFormat(float(np.sum(y))))
label = "{:s}".format(label)
ax0.step(x_all, y_all, label=label, lw=5, color="C0", where="mid")
ax0.plot(x, y_hat, linestyle="--", label="Fit", lw=5, color="C1")
ax0.plot([], [], ' ', label="$\\chi^{{2}}$/NDF = {:s}".format(LatexFormat(chi2/ndf)))
ax0.set_ylabel("Counts".format(xlabel, xlabel), fontsize=labelsize)
ax0.tick_params(axis="y", labelsize=ticksize)
if(y_limits[0] is None):
y_limits[0] = np.min(y[y>0])
ax0.set_ylim(y_limits)
ax0.set_xlim(x_limits)
if(title is not None):
ax0.set_title(title, fontsize=titlesize)
ax0.grid(which="both")
ax0.set_yscale("log")
ax0.legend(fontsize=legendsize)
ax1 = plt.subplot(gs[1], sharex = ax0)
ax1.scatter(x, (y - y_hat)/(y_hat_err), color="C1")
ax1.tick_params(axis="x", labelsize=ticksize, rotation=45)
ax1.tick_params(axis="y", labelsize=ticksize)
ax1.set_ylabel("Pull", fontsize=labelsize)
ax1.set_xlabel(xlabel, fontsize=labelsize)
ax1.set_ylim(residual_limits)
ax1.set_yticks(np.arange(int(np.floor(np.min(residual_limits))),
int(np.ceil(np.max(residual_limits)))))
ax1.axhline(y=-1.0, color="purple", lw=2.5, linestyle="--")
ax1.axhline(y=1.0, color="purple", lw=2.5, linestyle="--")
ax1.tick_params(axis="x", labelsize=ticksize)
ax1.tick_params(axis="y", labelsize=residualticksize)
ax1.grid(which="both")
fig.tight_layout()
if(not plot_in_bins):
mf = ScalarFormatter(useMathText=True)
mf.set_powerlimits((-2,2))
plt.gca().xaxis.set_major_formatter(mf)
offset = ax1.xaxis.get_offset_text()
offset.set_size(int(0.8*ticksize))
fig.subplots_adjust(hspace=.0)
if(save_directory is not None):
print("Saving figure to {0}...".format(save_directory))
fig.savefig(save_directory)
if(display):
plt.pause(0.01)
fig.show()
else:
plt.close(fig)
########################################################################################################
## Fit(prefit=False)
########################################################################################################
##
## VARIABLES:
## tau: slow time-constant of SiPM pulse in nanoseconds
## tau_err: optional error to allow tau to float
## t_0: pre-gate integration window in nanoseconds.
## t_0_err: optional error to allow t_0 to float
## tau_R: recovery time of SiPM in nanoseconds
## tau_R_err: optional error to allow tau_R to float
## t_gate: gate length in nanoseconds
## t_gate_err: optional error to allow t_gate to float
## bw: optional bin width if data is 1D
## truncate_nG: optional parameter to truncate distribution a certain fraction/number of estimated gain units before the estimated pedestal.
## peak_dist_nG: required minimum distance between peaks in gain units to register a peaks. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html
## peak_rel_height: Chooses the relative height at which the peak width is measured as a percentage of its prominence. 1.0 calculates the width of the peak at its lowest contour line while 0.5 evaluates at half the prominence height. Must be at least 0. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.peak_widths.html#scipy.signal.peak_widths
## bin_method: binning method for unbinned data. "knuth", "freedman" or "scott" rules are available.
## alpha_peaks: the maximum percentile of the distribution a found peak is allowed to have. This determines the threshold on the selected peak positions.
##
## Main fit function of PeakOTron, which performs the prefit and main fit using Minuit.
##
########################################################################################################
def __MinuitFit(self, f_DRM, A):
t_0 = time.time()
########################################################################################################
## Inititalise Minuit with hypothesis of no dark counts or afterpulsing.
########################################################################################################
m_DRM = Minuit(f_DRM,
Q_0 = self._prefit_values_bin["Q_0"],
G = self._prefit_values_bin["G"],
mu = self._prefit_values_bin["mu"],
lbda = self._prefit_values_bin["lbda"],
sigma_0 = self._prefit_values_bin["sigma_0"],
sigma_1 = self._prefit_values_bin["sigma_1"],
tau_Ap = self._prefit_values_bin["tau_Ap"],
p_Ap = epsilon(),
DCR = epsilon(),
tau = self._prefit_values_bin["tau"],
tau_R = self._prefit_values_bin["tau_R"],
t_gate = self._prefit_values_bin["t_gate"],
t_0 = self._prefit_values_bin["t_0"],
k_max = self._prefit_values_bin["k_max"],
k_max_dcr = self._prefit_values_bin["k_max_dcr"],
len_pad = self._prefit_values_bin["len_pad"],
A = A
)
m_DRM.errors["Q_0"] = self._prefit_errors_bin["Q_0"]
m_DRM.errors["G"] = self._prefit_errors_bin["G"]
m_DRM.errors["mu"] = self._prefit_errors_bin["mu"]
m_DRM.errors["lbda"] = self._prefit_errors_bin["lbda"]
m_DRM.errors["sigma_0"] = self._prefit_errors_bin["sigma_0"]
m_DRM.errors["sigma_1"] = self._prefit_errors_bin["sigma_1"]
m_DRM.errors["tau_Ap"] = self._prefit_errors_bin["tau_Ap"]
m_DRM.errors["p_Ap"] = self._prefit_errors_bin["p_Ap"]
m_DRM.errors["DCR"] = self._prefit_errors_bin["DCR"]
m_DRM.errors["A"] = np.sqrt(A)
########################################################################################################
## Set parameter limits
## These limits have been chosen to best constrain the possible values for each of the parameters.
## NOTE: we are working in bin widths, so the minimum allowed pedestal position is the first/second bin for Q_0, G, sigma_0, sigma_1
########################################################################################################
m_DRM.limits["G"] = (1, None)
m_DRM.limits["mu"] = (0.01, None)
m_DRM.limits["lbda"] = (1e-6, 1-1e-4)
m_DRM.limits["sigma_0"] = (0.1, None)
m_DRM.limits["sigma_1"] = (0.1, None)
m_DRM.limits["tau_Ap"] = (3, self._prefit_values_bin["t_gate"]/2)
m_DRM.limits["p_Ap"] = (1e-6, 1-1e-4)
m_DRM.limits["DCR"] = (1e-6, None)
m_DRM.limits["A"] = (self._prefit_values_bin["A"] - 3*np.sqrt(self._prefit_values_bin["A"]),
self._prefit_values_bin["A"] + 3*np.sqrt(self._prefit_values_bin["A"]))
# m_DRM.limits["A"] = (1,None)
########################################################################################################
## Allow input tau, t_gate, tau_R and t_0 to float if errors are known.
########################################################################################################
if(self._fit_kwargs["tau_err"] is None):
m_DRM.fixed["tau"]=True
else:
m_DRM.fixed["tau"]=False
m_DRM.errors["tau"] = self._prefit_errors_bin["tau"]
m_DRM.limits["tau"] = (max(epsilon(), self._prefit_values_bin["tau"]-self._prefit_errors_bin["tau"]),
max(epsilon(), self._prefit_values_bin["tau"]+self._prefit_errors_bin["tau"])
)
if(self._fit_kwargs["t_gate_err"] is None):
m_DRM.fixed["t_gate"]=True
else:
m_DRM.fixed["t_gate"]=False
m_DRM.errors["t_gate"] = self._prefit_errors_bin["t_gate"]
m_DRM.limits["t_gate"] = (max(epsilon(), self._prefit_values_bin["t_gate"]-self._prefit_errors_bin["t_gate"]),
max(epsilon(), self._prefit_values_bin["t_gate"]+self._prefit_errors_bin["t_gate"])
)
if(self._fit_kwargs["t_0_err"] is None):
m_DRM.fixed["t_0"]=True
else:
m_DRM.fixed["t_0"]=False
m_DRM.errors["t_0"] = self._prefit_errors_bin["t_0"]
m_DRM.limits["t_0"] = (max(epsilon(), self._prefit_values_bin["t_0"]-self._prefit_errors_bin["t_0"]),
max(epsilon(), self._prefit_values_bin["t_0"]+self._prefit_errors_bin["t_0"])
)
if(self._fit_kwargs["tau_R_err"] is None):
m_DRM.fixed["tau_R"]=True
else:
m_DRM.fixed["tau_R"]=False
m_DRM.errors["tau_R"] = self._prefit_errors_bin["tau_R"]
m_DRM.limits["tau_R"] = (max(epsilon(), self._prefit_values_bin["tau_R"]-self._prefit_errors_bin["tau_R"]),
max(epsilon(), self._prefit_values_bin["tau_R"]+self._prefit_errors_bin["tau_R"])
)
m_DRM.fixed["k_max"]=True
m_DRM.fixed["k_max_dcr"]=True
m_DRM.fixed["len_pad"]=True
########################################################################################################
##Perform fit under hypothesis there are no dark counts, afterpulses.
########################################################################################################
m_DRM.fixed["Q_0"] = False
m_DRM.fixed["G"] = False
m_DRM.fixed["mu"] = False
m_DRM.fixed["lbda"] = False
m_DRM.fixed["sigma_0"] = False
m_DRM.fixed["sigma_1"] = False
m_DRM.fixed["p_Ap"]=True
m_DRM.fixed["tau_Ap"]=True
m_DRM.fixed["DCR"]=True
m_DRM.fixed["A"]=False
if(self._verbose):
print("Fitting basic light...")
m_DRM.strategy=1
m_DRM.errordef=m_DRM.LIKELIHOOD
m_DRM.migrad(ncall = self._n_call_minuit,
iterate= self._n_iterations_minuit)
if(self._verbose):
print("First pass complete...")
print(m_DRM)
if(self._verbose):
print("Fixing basic light and fitting noise...")
m_DRM.fixed["Q_0"] = True
m_DRM.fixed["G"] = True
m_DRM.fixed["mu"] = True
m_DRM.fixed["lbda"] = True
m_DRM.fixed["sigma_0"] = True
m_DRM.fixed["sigma_1"] = True
m_DRM.fixed["p_Ap"]=False
m_DRM.fixed["tau_Ap"]=False
m_DRM.fixed["DCR"]=False
m_DRM.values["DCR"] = self._prefit_values_bin["DCR"]
m_DRM.values["p_Ap"] = self._prefit_values_bin["p_Ap"]
m_DRM.fixed["A"]=False
m_DRM.strategy=1
m_DRM.errordef=m_DRM.LIKELIHOOD
m_DRM.migrad(ncall = self._n_call_minuit,
iterate= self._n_iterations_minuit)
if(self._verbose):
print("Second pass complete...")
print(m_DRM)
if(self._verbose):
print("Fitting full model...")
m_DRM.fixed["Q_0"] = False
m_DRM.fixed["G"] = False
m_DRM.fixed["mu"] = False
m_DRM.fixed["lbda"] = False
m_DRM.fixed["sigma_0"] = False
m_DRM.fixed["sigma_1"] = False
m_DRM.fixed["p_Ap"]=False
m_DRM.fixed["tau_Ap"]=False
m_DRM.fixed["DCR"]=False
m_DRM.fixed["A"]=False
m_DRM.strategy=1
m_DRM.migrad(ncall = self._n_call_minuit,
iterate= self._n_iterations_minuit)
m_DRM.hesse()
t_1 = time.time()
time_fit = t_1 - t_0
if(self._verbose):
print("Prefit took {:3.3f} s, Fit took {:3.3f} s".format(self._time_prefit,
time_fit))
print(m_DRM)
return m_DRM, time_fit
def Fit(self, data, **kwargs):
########################################################################################################
## Preparation
########################################################################################################
##
## 1. Data dictionaries are initialised if previously filled.
## 2. tau, tauR, t_0 and t_gate are checked to ensure they are provided to the program.
## 3. Perform checks of input parameters to ensure no disallowed values
########################################################################################################
self.InitKwargs()
if not set(kwargs.keys()).issubset(self._fit_kwargs.keys()):
wrong_keys = list(set(kwargs.keys()).difference(set(self._fit_kwargs.keys())))
raise Exception ("Arguments {:s} not recognised.".format(",".join(wrong_keys)))
if not "tau" in kwargs.keys():
raise Exception ("Please provide 'tau' (slow component of SiPM pulse in nanoseconds).")
elif(not isinstance(kwargs["tau"], float) or kwargs["tau"]<=0):
raise Exception ("Please set 'tau' (slow component of SiPM pulse in nanoseconds) to a positive float.")
if not "tau_R" in kwargs.keys():
raise Exception ("Please provide 'tau_R' (recovery time of SiPM in nanoseconds).")
elif(not isinstance(kwargs["tau_R"], float) or kwargs["tau_R"]<=0):
raise Exception ("Please set 'tau' (slow component of SiPM pulse in nanoseconds) to a positive float.")
if not "t_0" in kwargs.keys():
raise Exception ("Please provide 't_0' (integration region before start of SiPM pulse integration gate in nanoseconds).")
elif(kwargs["t_0"]<=0):
raise Exception ("Please set 't_0' (integration region before start of SiPM pulse integration gate in nanoseconds) to a positive float.")
if not "t_gate" in kwargs.keys():
raise Exception ("Please provide 't_gate' (length of SiPM pulse integration gate) in nanoseconds.")
elif(kwargs["t_0"]<=0):
raise Exception ("Please set 't_0' (integration region before start of SiPM pulse integration gate in nanoseconds) to a positive float.")
self._fit_kwargs.update(kwargs)
if self._fit_kwargs["min_n_events_peak"] is not None:
if(not isinstance(self._fit_kwargs["min_n_events_peak"], int)):
raise Exception ("Please set max_n_events peak (max events to fit Gaus to peak) as an integer" )
else:
self._fit_kwargs["min_n_events_peak"] = self._min_n_events_peak
if self._fit_kwargs["truncate_nsigma0_up"] is not None:
if(not isinstance(self._fit_kwargs["truncate_nsigma0_up"], float)):
raise Exception ("Please set truncate_nsigma0_up (# of sigma0 before pedestal to for Chi2/NDF to be evaluated) as a positive float" )
elif(self._fit_kwargs["truncate_nsigma0_up"]<=0):
raise Exception ("Please set truncate_nsigma0_up (# of sigma0 before pedestal to for Chi2/NDF to be evaluated) as a positive float" )
if self._fit_kwargs["truncate_nsigma0_do"] is not None:
if(not isinstance(self._fit_kwargs["truncate_nsigma0_do"], float)):
raise Exception ("Please set truncate_nsigma0_do (max # of sigma0 before pedestal to for to be scanned) as a positive float between 0.0 and 4.0" )
elif(self._fit_kwargs["truncate_nsigma0_do"]<=0 or self._fit_kwargs["truncate_nsigma0_do"]>4):
raise Exception ("Please set truncate_nsigma0_do (# of sigma0 before pedestal to for Chi2/NDF to be evaluated) as a positive float between 0.0 and 4.0" )
if self._fit_kwargs["truncate_nsigma0_chi2scanlim"] is not None:
if(not isinstance(self._fit_kwargs["truncate_nsigma0_chi2scanlim"], float)):
raise Exception ("Please set truncate_nsigma0_chi2scanlim (minimum Chi2/NDF for a fit to be selected in a scan) as a positive float" )
elif(self._fit_kwargs["truncate_nsigma0_chi2scanlim"]<0):
raise Exception ("Please set truncate_nsigma0_chi2scanlim (minimum Chi2/NDF for a fit to be selected in a scan) as a positive float" )
if self._fit_kwargs["bw"] is not None:
if(not isinstance(self._fit_kwargs["bw"], float)):
raise Exception ("Please set bw (bin width) to a positive float." )
if(self._fit_kwargs["bw"]<=0):
raise Exception ("Please set bw (bin width) to a positive float." )
if self._fit_kwargs["bin_0"] is not None:
if(not isinstance(self._fit_kwargs["bin_0"], float)):
raise Exception ("Please set bin_0 (first bin position) to a float." )
if(self._verbose):
print("Prefitting...")
########################################################################################################
## Prefitting routine applied
########################################################################################################
t_0 = time.time()
self.GetPreFit(data, **self._fit_kwargs)
t_1 = time.time()
self._time_prefit = t_1-t_0
if(self._verbose):
print("Prefitting complete.")
_offset = self._hist_data["bc_min"]
_scale = self._hist_data["bw"]
if(self._fit_kwargs["prefit_only"]):
_prefit_values_temp = self._prefit_values_bin.copy()
_prefit_errors_temp = self._prefit_errors_bin.copy()
for key, value in _prefit_values_temp.items():
if(key=="Q_0"):
_prefit_values_temp[key] = _offset + value*_scale
elif(key=="G"):
_prefit_values_temp[key] = value*_scale
elif(key=="sigma_0"):
_prefit_values_temp[key] = value*_scale
elif(key=="sigma_1"):
_prefit_values_temp[key] = value*_scale
else:
_prefit_values_temp[key] = value
for key, value in _prefit_errors_temp.items():
if(key=="Q_0"):
_prefit_errors_temp[key] = value*_scale
elif(key=="G"):
_prefit_errors_temp[key] = value*_scale
elif(key=="sigma_0"):
_prefit_errors_temp[key] = value*_scale
elif(key=="sigma_1"):
_prefit_errors_temp[key] = value*_scale
else:
_prefit_errors_temp[key] = value
self._prefit_values.update(_prefit_values_temp)
self._prefit_errors.update(_prefit_errors_temp)
else:
########################################################################################################
## Truncate fit if supplied factor.
########################################################################################################
# lim_nsigma0 = self._prefit_values_bin["Q_0"] - self._fit_kwargs["truncate_nsigma0"]*self._prefit_values_bin["sigma_0"]
# if(self._verbose):
# print("Truncating fit, starting from {:3.3f} estimated sigma_0 from estimated pedestal ({:3.3f} bins) ...".format(self._fit_kwargs["truncate_nsigma0"], lim_nsigma0))
if(self._fit_kwargs["truncate_nsigma0_do"] is not None):
nsigma_0_scan_range = np.linspace(4.0,
self._fit_kwargs["truncate_nsigma0_do"],
np.floor((4-self._fit_kwargs["truncate_nsigma0_do"])/0.5).astype(int)+1)
if(self._verbose):
print("Performing full truncation scan. Evaluating nsigma0 = {0} before pedestal".format(nsigma_0_scan_range))
if(self._fit_kwargs["truncate_nsigma0_chi2scanlim"] is not None):
if(self._verbose):
print("Minimum Chi2/NDF for convergence: {:3.3f}".format(self._fit_kwargs["truncate_nsigma0_chi2scanlim"]))
else:
print("Best loss taken.")
if(self._fit_kwargs["truncate_nsigma0_up"] is not None):
_nsigma0hi = self._fit_kwargs["truncate_nsigma0_up"]
else:
_nsigma0hi = 2
if(self._verbose):
print("Evaluating up to nsigma0 = {:3.3f} after pedestal".format(_nsigma0hi))
m_DRMs = np.full(len(nsigma_0_scan_range), None)
chi2ndfs_range = np.full(len(nsigma_0_scan_range), None)
chi2ndfs_full = np.full(len(nsigma_0_scan_range), None)
time_fits = np.full(len(nsigma_0_scan_range), 0.0)
Q_0_prefit = self._prefit_values_bin["Q_0"]
sigma_0_prefit = self._prefit_values_bin["sigma_0"]
fit_iters = 0
for fit_i, _nsigma0low in enumerate(nsigma_0_scan_range):
fit_iters+=1
lim_nsigma0 = Q_0_prefit - _nsigma0low*sigma_0_prefit
print("Truncating fit, starting from {:3.3f} estimated sigma_0 from estimated pedestal ({:3.3f} bins) ...".format(_nsigma0low, lim_nsigma0))
cond_sel_nsigma0 = (self._hist_data["bin_numbers"]>lim_nsigma0)
f_DRM = BinnedLH(DRM,
self._hist_data["bin_numbers"][cond_sel_nsigma0],
self._hist_data["counts"][cond_sel_nsigma0], 1)
m_DRM, time_fit = self.__MinuitFit(f_DRM, np.sum(self._hist_data["counts"][cond_sel_nsigma0]))
chi2_range, ndof_range = self.GetChi2(pars = m_DRM.values.to_dict(),
n_sigma_lims = [_nsigma0low, _nsigma0hi]
)
red_chi2_range = chi2_range/ndof_range
chi2_full, ndof_full = self.GetChi2(pars = m_DRM.values.to_dict(),
n_sigma_lims = [_nsigma0low, None]
)
red_chi2_full = chi2_full/ndof_full
m_DRMs[fit_i] = m_DRM
chi2ndfs_range[fit_i] = red_chi2_range
chi2ndfs_full[fit_i] = red_chi2_full
time_fits[fit_i] = time_fit
if(self._fit_kwargs["truncate_nsigma0_chi2scanlim"] is not None):
if(red_chi2_range<self._fit_kwargs["truncate_nsigma0_chi2scanlim"]):
if(self._verbose):
print("Fit chi2/NDF less than limit. (Fit: {:3.3f} < Limit {:3.3f}. Scan complete".format(red_chi2,
self._fit_kwargs["truncate_nsigma0_chi2scanlim"]))
self._m_DRM = deepcopy(m_DRM)
self._time_fit = np.sum(time_fits)
self._hist_data["truncate_nsigma0_do_min"] = _nsigma0low
break
else:
if(self._verbose):
print("Fit chi2/NDF greater than limit. (Fit: {:3.3f} < Limit {:3.3f}. Performing next scan step...".format(red_chi2,
self._fit_kwargs["truncate_nsigma0_chi2scanlim"]))
if(fit_iters==len(nsigma_0_scan_range)):
fit_i_best = np.argmin(chi2ndfs_range)
if(self._verbose):
print("Scan Summary:")
for fit_i, _nsigma0 in enumerate(nsigma_0_scan_range):
if(fit_i==fit_i_best):
is_best = True
else:
is_best = False
print("n_sigma0 = {:3.3f}\tchi2/NDF in pedestal range = {:3.3f}\tBest = {:s}".format(_nsigma0, chi2ndfs_range[fit_i], str(is_best)))
print("Scan complete. Using best fit at {:3.3f}".format(nsigma_0_scan_range[fit_i_best]))
self._m_DRM = deepcopy(m_DRMs[fit_i_best])
self._time_fit = np.sum(time_fits)
#print(self._hist_data["truncate_nsigma0_do_min"] )
self._hist_data["truncate_nsigma0_do_min"] = nsigma_0_scan_range[fit_i_best]
#print(self._hist_data["truncate_nsigma0_do_min"] )
else:
f_DRM = BinnedLH(DRM,
self._hist_data["bin_numbers"],
self._hist_data["counts"], 1)
m_DRM, time_fit = self.__MinuitFit(f_DRM, np.sum(self._hist_data["counts"]))
self._time_fit = time_fit
self._m_DRM = deepcopy(m_DRM)
########################################################################################################
##Rescale fitted parameters back to their input units.
########################################################################################################
for key, value in self._m_DRM.values.to_dict().items():
if(key=="Q_0"):
self._fit_values[key] = _offset + value*_scale
elif(key=="G"):
self._fit_values[key] = value*_scale
elif(key=="sigma_0"):
self._fit_values[key] = value*_scale
elif(key=="sigma_1"):
self._fit_values[key] = value*_scale
else:
self._fit_values[key] = value
self._fit_values_bin[key] = value
for key, value in self._m_DRM.errors.to_dict().items():
if(key=="Q_0"):
self._fit_errors[key] = value*_scale
elif(key=="G"):
self._fit_errors[key] = value*_scale
elif(key=="sigma_0"):
self._fit_errors[key] = value*_scale
elif(key=="sigma_1"):
self._fit_errors[key] = value*_scale
else:
self._fit_errors[key] = value
self._fit_errors_bin[key] = value
_prefit_values_temp = self._prefit_values_bin.copy()
_prefit_errors_temp = self._prefit_errors_bin.copy()
for key, value in _prefit_values_temp.items():
if(key=="Q_0"):
_prefit_values_temp[key] = _offset + value*_scale
elif(key=="G"):
_prefit_values_temp[key] = value*_scale
elif(key=="sigma_0"):
_prefit_values_temp[key] = value*_scale
elif(key=="sigma_1"):
_prefit_values_temp[key] = value*_scale
else:
_prefit_values_temp[key] = value
for key, value in _prefit_errors_temp.items():
if(key=="Q_0"):
_prefit_errors_temp[key] = value*_scale
elif(key=="G"):
_prefit_errors_temp[key] = value*_scale
elif(key=="sigma_0"):
_prefit_errors_temp[key] = value*_scale
elif(key=="sigma_1"):
_prefit_errors_temp[key] = value*_scale
else:
_prefit_errors_temp[key] = value
self._prefit_values.update(_prefit_values_temp)
self._prefit_errors.update(_prefit_errors_temp)
########################################################################################################
## GetPreFit(prefit=False)
########################################################################################################
##
## VARIABLES:
## tau: slow time-constant of SiPM pulse in nanoseconds
## tau_err: optional error to allow tau to float
## t_0: pre-gate integration window in nanoseconds.
## t_0_err: optional error to allow t_0 to float
## tau_R: recovery time of SiPM in nanoseconds
## tau_R_err: optional error to allow tau_R to float
## t_gate: gate length in nanoseconds
## t_gate_err: optional error to allow t_gate to float
## bw: optional bin width if data is 1D
## truncate_nG: optional parameter to truncate distribution a certain fraction/number of estimated gain units before the estimated pedestal.
## peak_dist_nG: required minimum distance between peaks in gain units to register a peaks. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html
## peak_rel_height: Chooses the relative height at which the peak width is measured as a percentage of its prominence. 1.0 calculates the width of the peak at its lowest contour line while 0.5 evaluates at half the prominence height. Must be at least 0. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.peak_widths.html#scipy.signal.peak_widths
## bin_method: binning method for unbinned data. "knuth", "freedman" or "scott" rules are available.
## alpha_peaks: the maximum percentile of the distribution a found peak is allowed to have. This determines the threshold on the selected peak positions.
##
## Perform pre-fitting routine.
########################################################################################################
def GetPreFit(self, data, **kwargs):
########################################################################################################
##Get reconstructed data, bin width, number of bins and bin values
########################################################################################################
data, bw, nbins, bins = self.GetBins(data, kwargs["bw"], kwargs["bin_method"], kwargs["bin_0"])
counts, _ = np.histogram(data, bins=bins)
counts_error = np.where(counts<1, 1, np.sqrt(counts))
bin_centres = (bins[:-1] + bins[1:])/2.
nbins = len(bin_centres)
bin_numbers = np.arange(0, nbins)
self._len_pad = nbins
########################################################################################################
## Perform FFT and calculate Power Spectrum. The inverse of the frequency of the maximum amplitude yields the gain.
########################################################################################################
nbins_pad = nbins+2*nbins
counts_pad = np.pad(counts, (nbins, nbins), "constant", constant_values=0)
fhat = np.fft.fft(counts_pad) #computes the fft
psd = fhat * np.conj(fhat)/nbins_pad
fft_freq_orig = (1/nbins_pad) * np.arange(nbins_pad) #frequency array
idxs_half = np.arange(1, np.floor(nbins_pad/2), dtype=np.int32)
fft_amplitude = np.abs(psd[idxs_half])
fft_freq = fft_freq_orig[idxs_half]
try:
f_interp_fft = InterpolatedUnivariateSpline(fft_freq, fft_amplitude, k=4)
except:
raise Exception("Interpolation of FFT failed. Fit cannot proceed.")
try:
fft_cr_pts = f_interp_fft.derivative().roots()
except:
raise Exception("Root finding of FFT failed. Fit cannot proceed.")
if(len(fft_cr_pts)<1):
raise Exception("No roots found. Fit cannot proceed.")
fft_cr_vals = f_interp_fft(fft_cr_pts)
inv_G_fft = fft_cr_pts[np.argmax(fft_cr_vals)]
G_fft = 1/inv_G_fft
self._len_pad = np.around(G_fft, 0).astype(int)
########################################################################################################
#Estimate the position of the pedestal by minimising the GP gain w.r.t the FFT Gain
########################################################################################################
mu_unsub = HistMean(counts, bin_numbers)
sigma = HistStd(counts, bin_numbers)
gamma = HistSkew(counts, bin_numbers)
f_interp_counts = UnivariateSpline(bin_numbers, counts, w=1/counts_error, k=4)
try:
counts_cr_pts = f_interp_counts.derivative().roots()
except:
raise Exception("Root finding of spectrum failed. Fit cannot proceed.")
if(len(counts_cr_pts)<1):
raise Exception("No roots of spectrum found. Fit cannot proceed.")
counts_cr_vals = f_interp_counts(counts_cr_pts)
max_peak = counts_cr_pts[np.argmax(counts_cr_vals)]
f_pedestal = lambda Q_0: (GP_gain(mu_unsub-Q_0, sigma, gamma) - G_fft)**2
m_pedestal = Minuit(f_pedestal, Q_0 = max_peak)
m_pedestal.errordef=m_pedestal.LEAST_SQUARES
m_pedestal.limits["Q_0"] = (0, max_peak+epsilon())
m_pedestal.migrad()
Q_0_est = m_pedestal.values["Q_0"]
idx_peaks = np.sort(np.concatenate( [np.arange(max_peak-G_fft, -G_fft/2, -G_fft),
np.array([max_peak]),
np.arange(max_peak+G_fft, nbins+G_fft/2, G_fft),
] ))
idx_Q_0 = np.argmin(abs(idx_peaks-Q_0_est))
idx_peaks = idx_peaks[idx_Q_0:-1]
########################################################################################################
#Select maximum peak position, and step in units of gain up until either end of the spectrum
########################################################################################################
########################################################################################################
# - Select pedestal as closest peak position to estimated pedestal.
# - Select only peaks up to pedestal
# - Get midpoints and end values.
# - Linearly interpolate midpoints
# - Subtract 'background', and threshold to be a minimum of zero.
########################################################################################################
counts_bg = np.full(len(idx_peaks)-1, None)
idx_bg = np.full(len(idx_peaks)-1, None)
for i_peak in range(len(idx_peaks)-1):
last_peak= idx_peaks[i_peak]
next_peak= idx_peaks[i_peak+1]
mid_peak = idx_peaks[i_peak] + G_fft/2
cond_sel_interpeak = (counts_cr_pts>=last_peak) & (counts_cr_pts<=next_peak)
if(sum(cond_sel_interpeak)>0):
counts_cr_pts_interpeak = counts_cr_pts[cond_sel_interpeak]
counts_cr_values_interpeak = counts_cr_vals[cond_sel_interpeak]
idx_min_interpeak = np.argmin(counts_cr_values_interpeak)
count_min_interpeak = counts_cr_values_interpeak[idx_min_interpeak]
idx_min_interpeak = counts_cr_pts_interpeak[idx_min_interpeak]
else:
if(self._verbose):
print("No local minimum found in peak range. Using value at peak + gain/2")
idx_min_interpeak = mid_peak
count_min_interpeak = counts[np.argmin(abs(bin_numbers-mid_peak))]
idx_bg[i_peak] = idx_min_interpeak
counts_bg[i_peak] = count_min_interpeak
try:
arg_sort_bg = np.argsort(idx_bg)
idx_bg = idx_bg[arg_sort_bg]
counts_bg = counts_bg[arg_sort_bg]
first_nz_bin = bin_numbers[np.argmax(counts>0)]
counts_first_nz_bin = counts[np.argmax(counts>0)]
idx_bg = np.pad(idx_bg, (1,0), "constant", constant_values=first_nz_bin)
counts_bg = np.pad(counts_bg, (1,0), "constant", constant_values=counts_first_nz_bin)
f_interp_bg = interp1d(idx_bg, counts_bg, bounds_error=False, fill_value=(0,0))
counts_bg = f_interp_bg(bin_numbers).astype(float)
counts_bgsub = counts - counts_bg
counts_bgsub = np.where(counts_bgsub<0, 0, counts_bgsub)
counts_bgsub_error = np.where(counts_bgsub>1, np.sqrt(counts_bgsub), 1)
counts_bgsub_error = np.sqrt(counts_bgsub_error**2 + counts_error**2)
except:
if(self._verbose):
print("Background subtraction failed. Setting background subtracted spectrum to have no events.")
counts_bgsub = counts
counts_bg = np.zeros_like(counts)
counts_bgsub_error = counts_error
########################################################################################################
#Find maximum from spline fit of background subtracted
########################################################################################################
try:
f_interp_counts_bgsub = UnivariateSpline(bin_numbers,
counts_bgsub,
w = 1/counts_bgsub_error,
k=4)
except:
raise Exception("Interpolation of background-subtracted spectrum failed. Fit cannot proceed.")
try:
counts_cr_pts_bgsub = f_interp_counts_bgsub.derivative().roots()
except:
raise Exception("Root finding of background-subtracted spectrum failed. Fit cannot proceed.")
if(len(counts_cr_pts_bgsub)<1):
raise Exception("No roots found in background-subtracted spectrum. Fit cannot proceed.")
counts_cr_vals_bgsub = f_interp_counts_bgsub(counts_cr_pts_bgsub)
max_peak = counts_cr_pts_bgsub[np.argmax(counts_cr_vals_bgsub)]
idx_peaks = np.sort(np.concatenate( [np.arange(max_peak-G_fft, -G_fft/2, -G_fft),
np.array([max_peak]),
np.arange(max_peak+G_fft, nbins+G_fft/2, G_fft)] ))
########################################################################################################
# - Select pedestal as closest peak position to estimated pedestal.
# - Select only peaks up to pedestal
# - Update pedestal position.
########################################################################################################
idx_Q_0 = np.argmin(abs(idx_peaks-Q_0_est))
idx_peaks = idx_peaks[idx_Q_0:-1]
Q_0 = idx_peaks[0]
########################################################################################################
# - Find means and variances of all peaks in background-subtracted spectrum with greater than 100 events.
# - Fit line to means to get Gain/Pedestal - correction for binning
# - Fit line to variances to Get sigma0/sigma1
########################################################################################################
####################
#ITERATIVE GAUS FIT METHOD
#####################
idx_peak_mus = np.full(len(idx_peaks), np.nan)
idx_peak_dmus = np.full(len(idx_peaks), np.nan)
idx_peak_sigmas = np.full(len(idx_peaks), np.nan)
idx_peak_dsigmas = np.full(len(idx_peaks), np.nan)
idx_peak_vars = np.full(len(idx_peaks), np.nan)
idx_peak_dvars = np.full(len(idx_peaks), np.nan)
idx_peak_numbers= np.arange(0, len(idx_peaks)).astype(int)
f_gaus = lambda x, mu, sigma, A: A*norm.pdf(x, loc=mu, scale=sigma)
for p_i, peak in zip(idx_peak_numbers, idx_peaks):
min_rng_peak = peak - G_fft/2
max_rng_peak = peak + G_fft/2
cond_sel_peak = (bin_numbers>min_rng_peak) & (bin_numbers<max_rng_peak)
counts_bgsub_peak = counts_bgsub[cond_sel_peak]
bin_numbers_peak = bin_numbers[cond_sel_peak]
N_peak = np.around(np.sum(counts_bgsub_peak), 0).astype(int)
mu_peak = HistMean(counts_bgsub_peak, bin_numbers_peak)
sigma_peak = HistStd(counts_bgsub_peak, bin_numbers_peak)
dmu_peak = sigma_peak/np.sqrt(N_peak)
dsigma_peak = sigma_peak/np.sqrt(2*N_peak-2)
if(N_peak < kwargs["min_n_events_peak"]):
if(self._verbose):
print("Peak {:d} had {:d} events, less than {:d} event limit. Using mean and standard deviation...".format(p_i, N_peak, kwargs["min_n_events_peak"]))
else:
delta_mu = np.inf
delta_sigma = np.inf
try:
iter_count = 0
iter_failed = False
while(iter_count<10 and (delta_mu>0.01*G_fft and delta_sigma>0.01*G_fft)):
mu_mns2sigma = max(mu_peak-2*sigma_peak, min_rng_peak)
mu_pls2sigma = min(mu_peak+2*sigma_peak, max_rng_peak)
cond_sel_2sigma = (bin_numbers_peak>mu_mns2sigma) & (bin_numbers_peak<mu_pls2sigma)
if(sum(cond_sel_2sigma)<3):
iter_failed=True
break
f_gaus_peak = BinnedLH(f_gaus,
bin_numbers_peak[cond_sel_2sigma],
counts_bgsub_peak[cond_sel_2sigma ],
1
)
m_gaus_peak = Minuit(f_gaus_peak, mu=mu_peak, sigma=sigma_peak, A = sum(cond_sel_2sigma))
m_gaus_peak.limits['mu'] = (mu_mns2sigma, mu_pls2sigma)
m_gaus_peak.limits['sigma'] = (epsilon(), 2.0*sigma_peak)
m_gaus_peak.limits["A"] = (max(N_peak - 3*np.sqrt(N_peak), 1), N_peak + 3*np.sqrt(N_peak))
m_gaus_peak.errordef=m_pedestal.LIKELIHOOD
m_gaus_peak.migrad()
m_gaus_peak.hesse()
mu_peak_i = m_gaus_peak.values["mu"]
sigma_peak_i = m_gaus_peak.values["sigma"]
dmu_peak_i = m_gaus_peak.errors["mu"]
dsigma_peak_i = m_gaus_peak.errors["sigma"]
if(iter_count==0):
delta_mu = mu_peak_i
delta_sigma = sigma_peak_i
else:
delta_mu = abs(mu_peak - mu_peak_i)
delta_sigma = abs(sigma_peak - sigma_peak_i)
mu_peak = mu_peak_i
sigma_peak = sigma_peak_i
dmu_peak = dmu_peak_i
dsigma_peak = dsigma_peak_i
iter_count+=1
if(iter_failed):
if(self._verbose):
print("Iterative 2-sigma Gaus fit for Peak {:d} failed because the fit range went below 3 bins in size. Using mean and standard deviation.".format(p_i))
mu_peak = HistMean(counts_bgsub_peak, bin_numbers_peak)
sigma_peak = HistStd(counts_bgsub_peak, bin_numbers_peak)
dmu_peak = sigma_peak/np.sqrt(N_peak)
dsigma_peak = sigma_peak/np.sqrt(2*N_peak-2)
except:
if(self._verbose):
print("Iterative 2-sigma Gaus fit for Peak {:d} failed. Using mean and standard deviation...".format(p_i))
mu_peak = HistMean(counts_bgsub_peak, bin_numbers_peak)
sigma_peak = HistStd(counts_bgsub_peak, bin_numbers_peak)
dmu_peak = sigma_peak/np.sqrt(N_peak)
dsigma_peak = sigma_peak/np.sqrt(2*N_peak-2)
idx_peak_mus[p_i] = mu_peak
idx_peak_dmus[p_i] = dmu_peak
idx_peak_sigmas[p_i] = sigma_peak
idx_peak_dsigmas[p_i] = dsigma_peak
cond_sel_peaks = (~np.isnan(idx_peak_mus)) & (~np.isnan(idx_peak_sigmas)) & (~np.isnan(idx_peak_dmus)) & (~np.isnan(idx_peak_dsigmas)) & (idx_peak_sigmas>0.1)
if(sum(cond_sel_peaks)<3):
raise Exception("Less than three peaks found. Cannot proceed. Try reducing the min_n_events_peak variable.")
idx_peak_numbers = idx_peak_numbers[cond_sel_peaks]
idx_peak_mus = idx_peak_mus[cond_sel_peaks]
idx_peak_sigmas = idx_peak_sigmas[cond_sel_peaks]
idx_peak_dmus = idx_peak_dmus[cond_sel_peaks]
idx_peak_dsigmas = idx_peak_dsigmas[cond_sel_peaks]
idx_peak_vars = idx_peak_sigmas**2
idx_peak_dvars = 2*idx_peak_sigmas*idx_peak_dsigmas
f_mu_lin = HuberRegression(lambda x, G, Q_0: Linear(x, G, Q_0),
idx_peak_numbers,
idx_peak_mus,
idx_peak_dmus)
f_var_lin = HuberRegression(lambda x, var_1, var_0: Linear(x, var_1, var_0),
idx_peak_numbers,
idx_peak_vars,
idx_peak_dvars,
)
m_mu_lin = Minuit(f_mu_lin, Q_0 = Q_0, G = G_fft)
m_mu_lin.limits["Q_0"] = (epsilon(), nbins)
m_mu_lin.fixed["G"] = True
m_mu_lin.errordef = m_mu_lin.LEAST_SQUARES
m_mu_lin.migrad()
m_mu_lin.hesse()
m_var_lin = Minuit(f_var_lin, var_0 = 1, var_1 = 1)
m_var_lin.limits["var_0"] = (epsilon(), nbins)
m_var_lin.limits["var_1"] = (epsilon(), nbins)
m_var_lin.errordef = m_var_lin.LEAST_SQUARES
m_var_lin.migrad()
m_var_lin.hesse()
Q_0 = m_mu_lin.values["Q_0"]
dQ_0 = m_mu_lin.errors["Q_0"]
sigma_0 = np.sqrt(m_var_lin.values["var_0"])
dsigma_0 = (0.5/sigma_0)*m_var_lin.errors["var_0"]
sigma_1 = np.sqrt(m_var_lin.values["var_1"])
dsigma_1 = (0.5/sigma_0)*m_var_lin.errors["var_1"]
########################################################################################################
# Use final estimate of pedestal to determine the mu, lambda
# Calculate k_max by summing up the GP distribution probabilities (cumulative distribution function) for all peaks in
# the spectrum, and then finding the peak closest to the C_%
########################################################################################################
mu_GP = max(epsilon(), GP_mu(mu_unsub-Q_0, sigma, gamma))
lbda_GP = min(max(epsilon(), GP_lbda(mu_unsub-Q_0, sigma, gamma)), 1-epsilon())
k_max = len(idx_peaks)
if(self._max_n_peaks is not None):
if(k_max>self._max_n_peaks):
if(self._verbose):
print("Max # light peaks ({:d}) exceeded set maximum number of light peaks. Setting # peaks to {:d}".format(k_max, self._max_n_peaks))
k_max = self._max_n_peaks
elif(k_max<4):
if(self._verbose):
print("Max # light peaks ({:d}) less than 4, so cannot fit. Setting # peaks to 4".format(k_max))
k_max = 4
########################################################################################################
# - Convert to charge units: (K= Q - Q_0)/G
# - Estimate Probability density in a range of 0.1 about the minimum, K=0.5. N_DCR_min = mean(counts[0.45<K<0.55])/0.1
# - Get N_DCR = sum(counts[0<K<0.5]) + 0.5*sum(counts[0.5<counts<1])
# - calculate DCR
########################################################################################################
K = (bin_numbers - Q_0)/G_fft
dK = abs(K[1]-K[0])
################################
# In order to do the interpolated integral, we need to change co-ordinate
#bw = dx = 1
#dK = dK/dx = 1/G dx
#dx = G*dK
#so ∫counts.dx = G*∫counts.dK
#This is necessary for correct calculation of the error.
################################
cond_NN = (K>=0.45)&(K<=0.55)
cond_Nc = (K<0.5)
NN_sum = max(np.sum(counts[cond_NN]),1)
NN = np.mean(counts[cond_NN])
if(NN!=NN):
NN = float(counts[np.argmin(abs(K-0.5))])
NN_sum = float(NN)
Nc = max(float(np.sum(counts[cond_Nc])), 1)
R0p5 = NN/Nc
DCR = R0p5/(4*kwargs["tau"]*dK)
DCR = DCR*np.exp(DCR*kwargs["tau"])
dDCR = DCR*(1/np.sqrt(max(1, NN_sum)))
mu_DCR = DCR*(kwargs["t_gate"] + kwargs["t_0"])
if(DCR!=DCR or DCR<epsilon()):
if(self._verbose):
print("DCR returned invalid, or less than machine epsilon. Choosing DCR to be machine epsilon...")
DCR = epsilon()
dDCR = 1e-9
if(self._max_n_peaks_dcr is not None):
if(self._verbose):
print("Setting # dark peaks to {:d}".format( self._max_n_peaks_dcr))
k_max_dcr = self._max_n_peaks_dcr
elif(k_max_dcr<4):
if(self._verbose):
print("Max # dark peaks ({:d}) less than 4, so cannot fit. Setting # peaks to 4".format(k_max_dcr))
k_max_dcr = 4
self._hist_data["bw"] = bw
self._hist_data["bc_min"] = bin_centres[0]
self._hist_data["bins"]=bins
self._hist_data["nbins"]=nbins
self._hist_data["counts"] = counts
self._hist_data["counts_error"] = counts_error
self._hist_data["counts_bg"] = counts_bg
self._hist_data["counts_bgsub"] = counts_bgsub
self._hist_data["bin_centres"] = bin_centres
self._hist_data["bin_numbers"] = bin_numbers
self._hist_data["fft_amplitude"] = fft_amplitude
self._hist_data["fft_freq"] = fft_freq
self._hist_data["G_fft"] = G_fft
self._hist_data["Q_0_est"] = Q_0_est
self._hist_data["peaks"]["peak_position"] = idx_peaks
self._hist_data["peaks"]["peak_index"] = idx_peak_numbers
self._hist_data["peaks"]["peak_mean"] = idx_peak_mus
self._hist_data["peaks"]["peak_mean_error"] = idx_peak_dmus
self._hist_data["peaks"]["peak_variance"] = idx_peak_vars
self._hist_data["peaks"]["peak_variance_error"] = idx_peak_dvars
self._hist_data["peaks"]["peak_std_deviation"] = idx_peak_sigmas
self._hist_data["peaks"]["peak_std_deviation_error"] = idx_peak_dsigmas
self._prefit_values_bin["Q_0"] = Q_0
self._prefit_values_bin["G"] = G_fft
self._prefit_values_bin["mu"] = mu_GP
self._prefit_values_bin["lbda"] = lbda_GP
self._prefit_values_bin["sigma_0"] = sigma_0
self._prefit_values_bin["sigma_1"] = sigma_1
self._prefit_values_bin["DCR"] = DCR
self._prefit_values_bin["p_Ap"] = 0.1
self._prefit_values_bin["tau_Ap"] = 5.0
self._prefit_values_bin["tau"] = kwargs["tau"]
self._prefit_values_bin["t_0"] = kwargs["t_0"]
self._prefit_values_bin["t_gate"] = kwargs["t_gate"]
self._prefit_values_bin["tau_R"] = kwargs["tau_R"]
self._prefit_values_bin["k_max"]=k_max
self._prefit_values_bin["k_max_dcr"]=k_max_dcr
self._prefit_values_bin["len_pad"]=self._len_pad
self._prefit_values_bin["A"]=np.sum(counts)
self._prefit_errors_bin["Q_0"] = dQ_0
self._prefit_errors_bin["G"] = 1
self._prefit_errors_bin["mu"] = 0.1
self._prefit_errors_bin["lbda"] = 0.01
self._prefit_errors_bin["sigma_0"] = dsigma_0
self._prefit_errors_bin["sigma_1"] = dsigma_1
self._prefit_errors_bin["p_Ap"] = 0.05
self._prefit_errors_bin["tau_Ap"] = 0.05
self._prefit_errors_bin["tau"] = kwargs["tau_err"]
self._prefit_errors_bin["t_0"] = kwargs["t_0_err"]
self._prefit_errors_bin["tau_R"] = kwargs["tau_R_err"]
self._prefit_errors_bin["t_gate"] = kwargs["t_gate_err"]
self._prefit_errors_bin["DCR"] = dDCR
self._prefit_errors_bin["A"] = np.sqrt(np.sum(counts))