Skip to content
Snippets Groups Projects
Commit 357c0b31 authored by Jack Christopher Hutchinson Rolph's avatar Jack Christopher Hutchinson Rolph
Browse files

Update PeakOTron.py

parent 87a47487
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 93 additions and 49 deletions
PeakOTron.py
+ 93
− 49
View file @ 357c0b31
...@@ -284,6 +284,7 @@ class PeakOTron: ...@@ -284,6 +284,7 @@ class PeakOTron:
"bw":None, "bw":None,
"peak_dist_factor":0.8, "peak_dist_factor":0.8,
"peak_width_factor":0.5, "peak_width_factor":0.5,
"peak_ppf":0.99,
"prominence": 0.0, "prominence": 0.0,
"min_n_peak_evts": 100, "min_n_peak_evts": 100,
"kernel_kde":"gaussian", "kernel_kde":"gaussian",
...@@ -816,6 +817,7 @@ class PeakOTron: ...@@ -816,6 +817,7 @@ class PeakOTron:
ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="x", labelsize=self._plot_fontsize)
ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.set_yscale("log") ax.set_yscale("log")
ax.grid(which="both")
ax.legend(fontsize=max(10, self._plot_fontsize//1.5)) ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
...@@ -835,6 +837,7 @@ class PeakOTron: ...@@ -835,6 +837,7 @@ class PeakOTron:
ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="x", labelsize=self._plot_fontsize)
ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.set_yscale("log") ax.set_yscale("log")
ax.grid(which="both")
ax.set_ylim([1e-5, None]) ax.set_ylim([1e-5, None])
ax.legend(fontsize=max(10, self._plot_fontsize//1.5)) ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
...@@ -922,6 +925,7 @@ class PeakOTron: ...@@ -922,6 +925,7 @@ class PeakOTron:
ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="x", labelsize=self._plot_fontsize)
ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.set_yscale("log") ax.set_yscale("log")
ax.grid(which="both")
def PlotGenPois(self, ax): def PlotGenPois(self, ax):
...@@ -962,11 +966,14 @@ class PeakOTron: ...@@ -962,11 +966,14 @@ class PeakOTron:
ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="x", labelsize=self._plot_fontsize)
ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.legend(fontsize=max(10, self._plot_fontsize//1.5)) ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
ax.grid(which="both")
def PlotSigmaEst(self, ax): 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( 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"]),
self.LatexFormat(self._GD_data["fit_var0_err"]), self.LatexFormat(self._GD_data["fit_var0_err"]),
...@@ -975,15 +982,15 @@ class PeakOTron: ...@@ -975,15 +982,15 @@ class PeakOTron:
), fontsize=self._plot_fontsize) ), fontsize=self._plot_fontsize)
ax.errorbar(self._peak_data["n_peak_s"], ax.errorbar(self._peak_data["n_peak_s"][c],
self._peak_data["peak_var_s"], self._peak_data["peak_var_s"][c],
yerr=self._peak_data["peak_var_err_s"], yerr=self._peak_data["peak_var_err_s"][c],
label="Extracted Variances", label="Extracted Variances",
fmt="o") fmt="o")
ax.plot(self._peak_data["n_peak_s"], ax.plot(self._peak_data["n_peak_s"][c],
self.DummyLinear(self._peak_data["n_peak_s"], self.DummyLinear(self._peak_data["n_peak_s"][c],
self._GD_data["fit_var1"], self._GD_data["fit_var1"],
self._GD_data["fit_var0"] self._GD_data["fit_var0"]
), ),
...@@ -997,11 +1004,14 @@ class PeakOTron: ...@@ -997,11 +1004,14 @@ class PeakOTron:
ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="x", labelsize=self._plot_fontsize)
ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.legend(fontsize=max(10, self._plot_fontsize//1.5)) ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
ax.grid(which="both")
def PlotGainEst(self, ax): 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( ax.set_title("Estimated Gain \n $G$={0} $\pm$ {1} \n $PH_{{0}}$={2} $\pm$ {3}".format(
self.LatexFormat(self._GD_data["fit_G"]), self.LatexFormat(self._GD_data["fit_G"]),
self.LatexFormat(self._GD_data["fit_G_err"]), self.LatexFormat(self._GD_data["fit_G_err"]),
...@@ -1010,15 +1020,15 @@ class PeakOTron: ...@@ -1010,15 +1020,15 @@ class PeakOTron:
), fontsize=self._plot_fontsize) ), fontsize=self._plot_fontsize)
ax.errorbar(self._peak_data["n_peak_s"], ax.errorbar(self._peak_data["n_peak_s"][c],
self._peak_data["peak_mean_s"], self._peak_data["peak_mean_s"][c],
yerr=self._peak_data["peak_mean_err_s"], yerr=self._peak_data["peak_mean_err_s"][c],
label="Extracted Peak Means", label="Extracted Peak Means",
fmt="o") fmt="o")
ax.plot(self._peak_data["n_peak_s"], ax.plot(self._peak_data["n_peak_s"][c],
self.DummyLinear(self._peak_data["n_peak_s"], self.DummyLinear(self._peak_data["n_peak_s"][c],
self._GD_data["fit_G"], self._GD_data["fit_G"],
self._GD_data["fit_x_0"] self._GD_data["fit_x_0"]
), ),
...@@ -1032,6 +1042,7 @@ class PeakOTron: ...@@ -1032,6 +1042,7 @@ class PeakOTron:
ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="x", labelsize=self._plot_fontsize)
ax.tick_params(axis="y", labelsize=self._plot_fontsize) ax.tick_params(axis="y", labelsize=self._plot_fontsize)
ax.legend(fontsize=max(10, self._plot_fontsize//1.5)) ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
ax.grid(which="both")
...@@ -1112,7 +1123,7 @@ class PeakOTron: ...@@ -1112,7 +1123,7 @@ class PeakOTron:
fontsize=self._plot_fontsize) fontsize=self._plot_fontsize)
ax.set_xlabel("$PH_{norm, est}$", 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.set_ylabel("$\\frac{dp_{DCR}}{PH_{norm, est}}$", fontsize=self._plot_fontsize)
ax.tick_params(axis="x", labelsize=self._plot_fontsize) ax.tick_params(axis="x", labelsize=self._plot_fontsize)
...@@ -1187,7 +1198,7 @@ class PeakOTron: ...@@ -1187,7 +1198,7 @@ class PeakOTron:
if(title is not None): if(title is not None):
ax0.set_title(title, fontsize=titlesize) ax0.set_title(title, fontsize=titlesize)
ax0.grid(which="both")
ax0.legend(fontsize=legendsize) ax0.legend(fontsize=legendsize)
...@@ -1203,7 +1214,7 @@ class PeakOTron: ...@@ -1203,7 +1214,7 @@ class PeakOTron:
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=self._plot_fontsize) ax1.tick_params(axis="x", labelsize=self._plot_fontsize)
ax1.tick_params(axis="y", labelsize=self._plot_fontsize) ax1.tick_params(axis="y", labelsize=self._plot_fontsize)
ax1.grid(which="both")
fig.tight_layout() fig.tight_layout()
...@@ -1282,6 +1293,7 @@ class PeakOTron: ...@@ -1282,6 +1293,7 @@ class PeakOTron:
bw, bw,
peak_dist_factor, peak_dist_factor,
peak_width_factor, peak_width_factor,
peak_ppf,
pedestal, pedestal,
prominence, prominence,
bin_method, bin_method,
...@@ -1328,6 +1340,12 @@ class PeakOTron: ...@@ -1328,6 +1340,12 @@ class PeakOTron:
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 #Calculate histogram information
y_N, x = np.histogram(data, bins=bins, density=False) y_N, x = np.histogram(data, bins=bins, density=False)
...@@ -1342,6 +1360,7 @@ class PeakOTron: ...@@ -1342,6 +1360,7 @@ class PeakOTron:
y_N_s, x_s = np.histogram(data_s, bins=bins_s, density=False) y_N_s, x_s = np.histogram(data_s, bins=bins_s, density=False)
N_s = np.sum(y_N_s) N_s = np.sum(y_N_s)
...@@ -1354,6 +1373,7 @@ class PeakOTron: ...@@ -1354,6 +1373,7 @@ class PeakOTron:
try: try:
x_kde_s, y_kde_s, y_kde_err_s, _, = self.BootstrapKDE(data_s, x_kde_s, y_kde_s, y_kde_err_s, _, = self.BootstrapKDE(data_s,
...@@ -1437,12 +1457,21 @@ class PeakOTron: ...@@ -1437,12 +1457,21 @@ class PeakOTron:
self._peak_data["n_peak_s"] = [] self._peak_data["n_peak_s"] = []
self._peak_data["x_width_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 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)): 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)) 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): if(self._verbose):
print("Peak {:d} has {:d} events. Less than event threshold {:d}. Ignoring...".format(x_i, print("Peak {:d} has {:d} events. Less than event threshold {:d}. Ignoring...".format(x_i,
N_evt_peak, N_evt_peak,
...@@ -1646,7 +1675,7 @@ class PeakOTron: ...@@ -1646,7 +1675,7 @@ class PeakOTron:
self._peak_data["peak_var_s"]=[] self._peak_data["peak_var_s"]=[]
self._peak_data["peak_var_err_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"], 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["y_peak_s"],
self._peak_data["x_peak_s"], self._peak_data["x_peak_s"],
...@@ -1654,7 +1683,7 @@ class PeakOTron: ...@@ -1654,7 +1683,7 @@ class PeakOTron:
data_range = self.SelectRange(data_s, x_peak - x_width/2, x_peak + x_width/2) data_range = self.SelectRange(data_s, x_peak - x_width/2, x_peak + x_width/2)
try:
est_mean, std_err_mean, _ = self.Bootstrap(data_range, est_mean, std_err_mean, _ = self.Bootstrap(data_range,
np.mean, np.mean,
n_samples=n_bootstrap) n_samples=n_bootstrap)
...@@ -1667,18 +1696,32 @@ class PeakOTron: ...@@ -1667,18 +1696,32 @@ class PeakOTron:
self._peak_data["peak_mean_err_s"].append(std_err_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_s"].append(est_var)
self._peak_data["peak_var_err_s"].append(std_err_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_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_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_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"]) 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"] n_peaks_select = self._peak_data["n_peak_s"][cond_peak]
peaks_var_select = self._peak_data["peak_var_s"] peaks_var_select = self._peak_data["peak_var_s"][cond_peak]
peaks_var_err_select = self._peak_data["peak_var_err_s"] peaks_var_err_select = self._peak_data["peak_var_err_s"][cond_peak]
peaks_mean_select = self._peak_data["peak_mean_s"] peaks_mean_select = self._peak_data["peak_mean_s"][cond_peak]
peaks_mean_err_select = self._peak_data["peak_mean_err_s"] peaks_mean_err_select = self._peak_data["peak_mean_err_s"][cond_peak]
fit_lin_var = HuberRegression(self.DummyLinear, fit_lin_var = HuberRegression(self.DummyLinear,
n_peaks_select, n_peaks_select,
...@@ -1957,6 +2000,7 @@ class PeakOTron: ...@@ -1957,6 +2000,7 @@ class PeakOTron:
bw, bw,
peak_dist_factor, peak_dist_factor,
peak_width_factor, peak_width_factor,
peak_ppf,
prominence, prominence,
min_n_peak_evts, min_n_peak_evts,
kernel_kde, kernel_kde,
...@@ -1977,11 +2021,11 @@ class PeakOTron: ...@@ -1977,11 +2021,11 @@ class PeakOTron:
if(self._verbose): if(self._verbose):
print("Getting histogram and peak values...") print("Getting histogram and peak values...")
self.GetHistAndPeaks( self.GetHistAndPeaks(data,
data,
bw, bw,
peak_dist_factor, peak_dist_factor,
peak_width_factor, peak_width_factor,
peak_ppf,
pedestal, pedestal,
prominence, prominence,
bin_method, bin_method,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment