diff --git a/simulations/11_pulsed_timing.py b/simulations/11_pulsed_timing.py index 98ba4b8..1a91bba 100755 --- a/simulations/11_pulsed_timing.py +++ b/simulations/11_pulsed_timing.py @@ -5,7 +5,7 @@ from lib import util from scipy import signal, interpolate, stats import matplotlib.pyplot as plt import numpy as np -from itertools import zip_longest +from itertools import zip_longest, pairwise import h5py from copy import deepcopy @@ -232,6 +232,7 @@ if __name__ == "__main__": h5_cache_fname = f'11_pulsed_timing.hdf5' time_accuracies = np.zeros(len(snr_factors)) + mask_counts = np.zeros(len(snr_factors)) for k, snr_sigma_factor in tqdm(enumerate(snr_factors)): # Read in cached time residuals if True: @@ -448,65 +449,90 @@ if __name__ == "__main__": # Make a plot of the time residuals if N_residuals > 1: - time_accuracies[k] = np.std(time_residuals[:N_residuals]) + time_residuals = time_residuals[:N_residuals] - hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step') - fig, ax = plt.subplots() - ax.set_title( - "Template Correlation Lag finding" - + f"\n template dt: {template_dt*1e3: .1e}ps" - + f"; antenna dt: {antenna_dt: .1e}ns" - + f"; noise_factor: {noise_sigma_factor: .1e}" - ) - ax.set_xlabel("Time Residual [ns]") - ax.set_ylabel("#") + for i in range(1 + cut_wrong_peak_matches): + mask_count = 0 - counts, bins, _patches = ax.hist(time_residuals, **hist_kwargs) - if True: # fit gaussian to histogram - min_x = min(time_residuals) - max_x = max(time_residuals) + if i==1: # if cut_wrong_peak_matches: + wrong_peak_condition = lambda t_res: abs(t_res) > antenna_dt*4 - dx = bins[1] - bins[0] - scale = len(time_residuals) * dx + mask = wrong_peak_condition(time_residuals) - xs = np.linspace(min_x, max_x) + mask_count = np.count_nonzero(mask) - # do the fit - name = "Norm" - param_names = [ "$\\mu$", "$\\sigma$" ] - distr_func = stats.norm + print("Masking {} residuals".format(mask_count)) + time_residuals = time_residuals[~mask] - label = name +"(" + ','.join(param_names) + ')' + if not mask_count: + print("Continuing") + continue - # plot - fit_params = distr_func.fit(time_residuals) - fit_ys = scale * distr_func.pdf(xs, *fit_params) - ax.plot(xs, fit_ys, label=label) + time_accuracies[k] = np.std(time_residuals) + mask_counts[k] = mask_count + + hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step') + fig, ax = plt.subplots() + ax.set_title( + "Template Correlation Lag finding" + + f"\n template dt: {template_dt: .1e}ns" + + f"; antenna dt: {antenna_dt: .1e}ns" + + ";" if not mask_count else "\n" + + f"snr_factor: {snr_sigma_factor: .1e}" + + "" if not mask_count else f"; N_masked: {mask_count}" + ) + ax.set_xlabel("Time Residual [ns]") + ax.set_ylabel("#") + + counts, bins, _patches = ax.hist(time_residuals, **hist_kwargs) + if True: # fit gaussian to histogram + min_x = min(time_residuals) + max_x = max(time_residuals) + + dx = bins[1] - bins[0] + scale = len(time_residuals) * dx + + xs = np.linspace(min_x, max_x) + + # do the fit + name = "Norm" + param_names = [ "$\\mu$", "$\\sigma$" ] + distr_func = stats.norm + + label = name +"(" + ','.join(param_names) + ')' + + # plot + fit_params = distr_func.fit(time_residuals) + fit_ys = scale * distr_func.pdf(xs, *fit_params) + ax.plot(xs, fit_ys, label=label) + + # chisq + ct = np.diff(distr_func.cdf(bins, *fit_params))*np.sum(counts) + if True: + ct *= np.sum(counts)/np.sum(ct) + c2t = stats.chisquare(counts, ct, ddof=len(fit_params)) + chisq_strs = [ + f"$\\chi^2$/dof = {c2t[0]: .2g}/{len(fit_params)}" + ] + + # text on plot + text_str = "\n".join( + [label] + + + [ f"{param} = {value: .2e}" for param, value in zip_longest(param_names, fit_params, fillvalue='?') ] + + + chisq_strs + ) + + ax.text( *(0.02, 0.95), text_str, fontsize=12, ha='left', va='top', transform=ax.transAxes) + + if mask_count: + fig.savefig(f"figures/11_time_residual_hist_tdt{template_dt:0.1e}_n{snr_sigma_factor:.1e}_masked.pdf") + else: + fig.savefig(f"figures/11_time_residual_hist_tdt{template_dt:0.1e}_n{snr_sigma_factor:.1e}.pdf") - # chisq - ct = np.diff(distr_func.cdf(bins, *fit_params))*np.sum(counts) if True: - ct *= np.sum(counts)/np.sum(ct) - c2t = stats.chisquare(counts, ct, ddof=len(fit_params)) - chisq_strs = [ - f"$\\chi^2$/dof = {c2t[0]: .2g}/{len(fit_params)}" - ] - - # text on plot - text_str = "\n".join( - [label] - + - [ f"{param} = {value: .2e}" for param, value in zip_longest(param_names, fit_params, fillvalue='?') ] - + - chisq_strs - ) - - ax.text( *(0.02, 0.95), text_str, fontsize=12, ha='left', va='top', transform=ax.transAxes) - - fig.savefig(f"figures/11_time_residual_hist_tdt{template_dt:0.1e}_n{noise_sigma_factor: .1e}.pdf") - - if True: - plt.close(fig) + plt.close(fig) # SNR time accuracy plot if True: @@ -526,7 +552,17 @@ if __name__ == "__main__": ax.set_yscale('log') # plot the values - ax.plot(np.asarray(snr_factors), time_accuracies, ls='none', marker='o') + l = None + for j, mask_threshold in enumerate(pairwise([np.inf, 250, 50, 1, 0])): + kwargs = dict( + ls='none', + marker=['^', 'v','8', 'o',][j], + color=None if l is None else l[0].get_color(), + ) + mask = mask_counts >= mask_threshold[1] + mask &= mask_counts < mask_threshold[0] + + l = ax.plot(snr_factors[mask], time_accuracies[mask], **kwargs) if True: # limit y-axis to 1e0 ax.set_ylim([None, 1e1])