diff --git a/PeakOTron.py b/PeakOTron.py index b58fce335ed1e554fce262007ca2b8999307f20a..773ff6f150ce9b399d576edcd2792c0e2f73faf2 100644 --- a/PeakOTron.py +++ b/PeakOTron.py @@ -284,6 +284,7 @@ class PeakOTron: "bw":None, "peak_dist_factor":0.8, "peak_width_factor":0.5, + "peak_ppf":0.99, "prominence": 0.0, "min_n_peak_evts": 100, "kernel_kde":"gaussian", @@ -625,7 +626,7 @@ class PeakOTron: def ApplyDCR(self, x, y_light, x_0, G, DCR, tgate, t_0, tau, r_fast, lbda, dx, k_dcr_low, k_dcr_hi): mu_dcr = DCR*(abs(t_0) + tgate) - + x_lin = np.arange(0, len(x)) x_pad_lin = np.arange(-self._len_DCR_pad, len(x)+self._len_DCR_pad) f_lin = interp1d(x, x_lin, fill_value='extrapolate') @@ -816,6 +817,7 @@ class PeakOTron: ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.set_yscale("log") + ax.grid(which="both") ax.legend(fontsize=max(10, self._plot_fontsize//1.5)) @@ -835,6 +837,7 @@ class PeakOTron: ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.set_yscale("log") + ax.grid(which="both") ax.set_ylim([1e-5, None]) ax.legend(fontsize=max(10, self._plot_fontsize//1.5)) @@ -922,6 +925,7 @@ class PeakOTron: ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.set_yscale("log") + ax.grid(which="both") def PlotGenPois(self, ax): @@ -962,11 +966,14 @@ class PeakOTron: ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.legend(fontsize=max(10, self._plot_fontsize//1.5)) + ax.grid(which="both") def PlotSigmaEst(self, ax): + c = self._peak_data["cond_valid_peak"] + ax.set_title("Estimated Peak Variance \n $\\sigma^{{2}}_{{0}}$={0} $\pm$ {1} \n $\\sigma^{{2}}_{{1}}$={2} $\pm$ {3}".format( self.LatexFormat(self._GD_data["fit_var0"]), self.LatexFormat(self._GD_data["fit_var0_err"]), @@ -975,15 +982,15 @@ class PeakOTron: ), fontsize=self._plot_fontsize) - ax.errorbar(self._peak_data["n_peak_s"], - self._peak_data["peak_var_s"], - yerr=self._peak_data["peak_var_err_s"], + ax.errorbar(self._peak_data["n_peak_s"][c], + self._peak_data["peak_var_s"][c], + yerr=self._peak_data["peak_var_err_s"][c], label="Extracted Variances", fmt="o") - ax.plot(self._peak_data["n_peak_s"], - self.DummyLinear(self._peak_data["n_peak_s"], + ax.plot(self._peak_data["n_peak_s"][c], + self.DummyLinear(self._peak_data["n_peak_s"][c], self._GD_data["fit_var1"], self._GD_data["fit_var0"] ), @@ -997,10 +1004,13 @@ class PeakOTron: ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.legend(fontsize=max(10, self._plot_fontsize//1.5)) + ax.grid(which="both") def PlotGainEst(self, ax): + + c = self._peak_data["cond_valid_peak"] ax.set_title("Estimated Gain \n $G$={0} $\pm$ {1} \n $PH_{{0}}$={2} $\pm$ {3}".format( self.LatexFormat(self._GD_data["fit_G"]), @@ -1010,15 +1020,15 @@ class PeakOTron: ), fontsize=self._plot_fontsize) - ax.errorbar(self._peak_data["n_peak_s"], - self._peak_data["peak_mean_s"], - yerr=self._peak_data["peak_mean_err_s"], + ax.errorbar(self._peak_data["n_peak_s"][c], + self._peak_data["peak_mean_s"][c], + yerr=self._peak_data["peak_mean_err_s"][c], label="Extracted Peak Means", fmt="o") - ax.plot(self._peak_data["n_peak_s"], - self.DummyLinear(self._peak_data["n_peak_s"], + ax.plot(self._peak_data["n_peak_s"][c], + self.DummyLinear(self._peak_data["n_peak_s"][c], self._GD_data["fit_G"], self._GD_data["fit_x_0"] ), @@ -1032,6 +1042,7 @@ class PeakOTron: ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.legend(fontsize=max(10, self._plot_fontsize//1.5)) + ax.grid(which="both") @@ -1112,7 +1123,7 @@ class PeakOTron: fontsize=self._plot_fontsize) ax.set_xlabel("$PH_{norm, est}$", fontsize=self._plot_fontsize) - + ax.grid(which="both") ax.set_ylabel("$\\frac{dp_{DCR}}{PH_{norm, est}}$", fontsize=self._plot_fontsize) ax.tick_params(axis="x", labelsize=self._plot_fontsize) @@ -1187,7 +1198,7 @@ class PeakOTron: if(title is not None): ax0.set_title(title, fontsize=titlesize) - + ax0.grid(which="both") ax0.legend(fontsize=legendsize) @@ -1203,7 +1214,7 @@ class PeakOTron: ax1.axhline(y=1.0, color="purple", lw=2.5, linestyle="--") ax1.tick_params(axis="x", labelsize=self._plot_fontsize) ax1.tick_params(axis="y", labelsize=self._plot_fontsize) - + ax1.grid(which="both") fig.tight_layout() @@ -1282,6 +1293,7 @@ class PeakOTron: bw, peak_dist_factor, peak_width_factor, + peak_ppf, pedestal, prominence, bin_method, @@ -1327,6 +1339,12 @@ class PeakOTron: bins_s = data_s.min() + bw_s * np.arange(nbins_s + 1) + + + if(nbins<50): + print("Binning with bw = {:3.3E} produced only {:d} bins, less than limit of 50 bins. Setting fit fail status and continuing...".format(bw, nbins)) + self._failed = True + return #Calculate histogram information @@ -1341,6 +1359,7 @@ class PeakOTron: x = (x[:-1] + x[1:])/2. + y_N_s, x_s = np.histogram(data_s, bins=bins_s, density=False) @@ -1351,6 +1370,7 @@ class PeakOTron: y_err_s = y_s*np.sqrt(1/y_N_s + 1/N_s) x_s = (x_s[:-1] + x_s[1:])/2. + @@ -1437,12 +1457,21 @@ class PeakOTron: self._peak_data["n_peak_s"] = [] self._peak_data["x_width_s"] = [] + + if(peak_ppf is not None): + peak_threshold = self.GPoisPPF(peak_ppf, est_mu, est_lbda, self._max_n_peaks) + if(self._verbose): + print("For {:3.3}% of distribution, accepting all peaks below {:d}".format(peak_ppf*100.0, peak_threshold)) + else: + peak_threshold = 0 + print("Accepting all peaks.") + n_valid_peaks = 0 for x_i, (y_peak, y_peak_err, x_peak, x_width) in enumerate(zip(y_peaks_s, y_peaks_err_s, x_peaks_s, x_widths_s)): N_evt_peak = len(self.SelectRange(data_s, x_peak-x_width/2, x_peak+x_width/2)) - if(N_evt_peak<min_n_peak_evts): + if(x_i>peak_threshold and N_evt_peak<min_n_peak_evts): if(self._verbose): print("Peak {:d} has {:d} events. Less than event threshold {:d}. Ignoring...".format(x_i, N_evt_peak, @@ -1646,7 +1675,7 @@ class PeakOTron: self._peak_data["peak_var_s"]=[] self._peak_data["peak_var_err_s"]=[] - + cond_peak = np.ones(len(peak_val), bool) for x_i, y_peak, x_peak, x_width in zip(self._peak_data["n_peak_s"], self._peak_data["y_peak_s"], self._peak_data["x_peak_s"], @@ -1654,31 +1683,45 @@ class PeakOTron: data_range = self.SelectRange(data_s, x_peak - x_width/2, x_peak + x_width/2) - - est_mean, std_err_mean, _ = self.Bootstrap(data_range, - np.mean, - n_samples=n_bootstrap) - - est_var, std_err_var, _ = self.Bootstrap(data_range, - np.var, - n_samples=n_bootstrap) - - self._peak_data["peak_mean_s"].append(est_mean) - self._peak_data["peak_mean_err_s"].append(std_err_mean) - self._peak_data["peak_var_s"].append(est_var) - self._peak_data["peak_var_err_s"].append(std_err_var) - - + try: + est_mean, std_err_mean, _ = self.Bootstrap(data_range, + np.mean, + n_samples=n_bootstrap) + + est_var, std_err_var, _ = self.Bootstrap(data_range, + np.var, + n_samples=n_bootstrap) + + self._peak_data["peak_mean_s"].append(est_mean) + self._peak_data["peak_mean_err_s"].append(std_err_mean) + self._peak_data["peak_var_s"].append(est_var) + self._peak_data["peak_var_err_s"].append(std_err_var) + cond_peak[x_i]=True + except: + cond_peak[x_i]=False + + + self._peak_data["peak_mean_s"].append(None) + self._peak_data["peak_mean_err_s"].append(None) + self._peak_data["peak_var_s"].append(None) + self._peak_data["peak_var_err_s"].append(None) + + if(sum(cond_peak)<3): + print("Fewer than 3 peaks had adequate statistics to determine peak width. Setting fit fail status...") + self._failed=True + return + + self._peak_data["cond_valid_peak"] = cond_peak self._peak_data["peak_mean_s"] = np.asarray(self._peak_data["peak_mean_s"]) self._peak_data["peak_mean_err_s"] = np.asarray(self._peak_data["peak_mean_err_s"]) self._peak_data["peak_var_s"] = np.asarray(self._peak_data["peak_var_s"]) self._peak_data["peak_var_err_s"] = np.asarray(self._peak_data["peak_var_err_s"]) - n_peaks_select = self._peak_data["n_peak_s"] - peaks_var_select = self._peak_data["peak_var_s"] - peaks_var_err_select = self._peak_data["peak_var_err_s"] - peaks_mean_select = self._peak_data["peak_mean_s"] - peaks_mean_err_select = self._peak_data["peak_mean_err_s"] + n_peaks_select = self._peak_data["n_peak_s"][cond_peak] + peaks_var_select = self._peak_data["peak_var_s"][cond_peak] + peaks_var_err_select = self._peak_data["peak_var_err_s"][cond_peak] + peaks_mean_select = self._peak_data["peak_mean_s"][cond_peak] + peaks_mean_err_select = self._peak_data["peak_mean_err_s"][cond_peak] fit_lin_var = HuberRegression(self.DummyLinear, n_peaks_select, @@ -1957,6 +2000,7 @@ class PeakOTron: bw, peak_dist_factor, peak_width_factor, + peak_ppf, prominence, min_n_peak_evts, kernel_kde, @@ -1977,19 +2021,19 @@ class PeakOTron: if(self._verbose): print("Getting histogram and peak values...") - self.GetHistAndPeaks( - data, - bw, - peak_dist_factor, - peak_width_factor, - pedestal, - prominence, - bin_method, - n_bootstrap, - n_bootstrap_kde, - bandwidth_kde, - kernel_kde, - min_n_peak_evts) + self.GetHistAndPeaks(data, + bw, + peak_dist_factor, + peak_width_factor, + peak_ppf, + pedestal, + prominence, + bin_method, + n_bootstrap, + n_bootstrap_kde, + bandwidth_kde, + kernel_kde, + min_n_peak_evts) if(not self._failed): if(self._verbose):