From 59feab014e07df20dbbeae83d8d3397097ac9525 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Wed, 26 Apr 2023 17:04:34 +0200 Subject: [PATCH] Pulse: snr plot multiple template_dt curves --- simulations/11_pulsed_timing.py | 217 +++++++++++++++++--------------- 1 file changed, 114 insertions(+), 103 deletions(-) diff --git a/simulations/11_pulsed_timing.py b/simulations/11_pulsed_timing.py index a601abc..9da6cb4 100755 --- a/simulations/11_pulsed_timing.py +++ b/simulations/11_pulsed_timing.py @@ -401,11 +401,14 @@ if __name__ == "__main__": matplotlib.use('Agg') bp_freq = (30e-3, 80e-3) # GHz - template_dt = 5e-2 # ns interp_template_dt = 5e-5 # ns template_length = 200 # ns + antenna_dt = 2 # ns + antenna_timelength = 1024 # ns + N_residuals = 50*3 if len(sys.argv) < 2 else int(sys.argv[1]) + template_dts = np.array([antenna_dt, 5e-1, 5e-2]) # ns snr_factors = np.concatenate( # 1/noise_amplitude factor ( #[0.25, 0.5, 0.75], @@ -415,8 +418,6 @@ if __name__ == "__main__": ), axis=None, dtype=float) - antenna_dt = 2 # ns - antenna_timelength = 1024 # ns cut_wrong_peak_matches = True normalise_noise = False @@ -454,112 +455,113 @@ if __name__ == "__main__": if True: plt.close(fig) - # - # Create the template - # This is sampled at a lower samplerate than the interpolation template - # - template, _ = create_template(dt=template_dt, timelength=template_length, bp_freq=bp_freq, name='Template') # # Find time accuracies as a function of signal strength # - time_accuracies = np.zeros(len(snr_factors)) - mask_counts = np.zeros(len(snr_factors)) - for k, snr_sigma_factor in tqdm(enumerate(snr_factors)): + time_accuracies = np.zeros((len(template_dts), len(snr_factors))) + mask_counts = np.zeros_like(time_accuracies) + for l, template_dt in tqdm(enumerate(template_dts)): - time_residuals = get_time_residuals_for_template( - N_residuals, template, interpolation_template=interp_template, - antenna_dt=antenna_dt, antenna_timelength=antenna_timelength, - snr_sigma_factor=snr_sigma_factor, bp_freq=bp_freq, normalise_noise=normalise_noise, - h5_cache_fname=h5_cache_fname, rng=rng, tqdm=tqdm) + # Create the template + # This is sampled at a lower samplerate than the interpolation template + template, _ = create_template(dt=template_dt, timelength=template_length, bp_freq=bp_freq, name='Template') - print()# separating tqdm - print()# separating tqdm + for k, snr_sigma_factor in tqdm(enumerate(snr_factors)): - # Make a plot of the time residuals - if N_residuals > 1: - for i in range(1 + cut_wrong_peak_matches): - mask_count = 0 + # get the time residuals + time_residuals = get_time_residuals_for_template( + N_residuals, template, interpolation_template=interp_template, + antenna_dt=antenna_dt, antenna_timelength=antenna_timelength, + snr_sigma_factor=snr_sigma_factor, bp_freq=bp_freq, normalise_noise=normalise_noise, + h5_cache_fname=h5_cache_fname, rng=rng, tqdm=tqdm) - if i==1: # if cut_wrong_peak_matches: - wrong_peak_condition = lambda t_res: abs(t_res) > antenna_dt*4 + print()# separating tqdm + print()# separating tqdm - mask = wrong_peak_condition(time_residuals) + # Make a plot of the time residuals + if N_residuals > 1: + for i in range(1 + cut_wrong_peak_matches): + mask_count = 0 - mask_count = np.count_nonzero(mask) + if i==1: # if cut_wrong_peak_matches: + wrong_peak_condition = lambda t_res: abs(t_res) > antenna_dt*4 - print("Masking {} residuals".format(mask_count)) - time_residuals = time_residuals[~mask] + mask = wrong_peak_condition(time_residuals) - if not mask_count: - print("Continuing") - continue + mask_count = np.count_nonzero(mask) - time_accuracies[k] = np.std(time_residuals) - mask_counts[k] = mask_count + print("Masking {} residuals".format(mask_count)) + time_residuals = time_residuals[~mask] - 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("#") + if not mask_count: + continue - 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) + time_accuracies[l, k] = np.std(time_residuals) + mask_counts[l, k] = mask_count - dx = bins[1] - bins[0] - scale = len(time_residuals) * dx + 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("#") - xs = np.linspace(min_x, max_x) + 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) - # do the fit - name = "Norm" - param_names = [ "$\\mu$", "$\\sigma$" ] - distr_func = stats.norm + dx = bins[1] - bins[0] + scale = len(time_residuals) * dx - label = name +"(" + ','.join(param_names) + ')' + xs = np.linspace(min_x, max_x) - # plot - fit_params = distr_func.fit(time_residuals) - fit_ys = scale * distr_func.pdf(xs, *fit_params) - ax.plot(xs, fit_ys, label=label) + # 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) - - 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") - - if True: - plt.close(fig) + plt.close(fig) # SNR time accuracy plot if True: @@ -567,10 +569,11 @@ if __name__ == "__main__": ax.set_title(f"Template matching SNR vs time accuracy") ax.set_xlabel("Signal to Noise Factor") ax.set_ylabel("Time Accuracy [ns]") + ax.grid() ax.legend(title="\n".join([ f"N={N_residuals}", - f"template_dt={template_dt:0.1e}ns", + #f"template_dt={template_dt:0.1e}ns", f"antenna_dt={antenna_dt:0.1e}ns", ])) @@ -578,28 +581,36 @@ if __name__ == "__main__": ax.set_xscale('log') ax.set_yscale('log') - # plot the values - 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] + # plot the values per template_dt slice + template_dt_colors = [None]*len(template_dts) + for k, template_dt in enumerate(template_dts): - l = ax.plot(snr_factors[mask], time_accuracies[mask], **kwargs) + # indicate masking values + 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 template_dt_colors[k] is None else template_dt_colors[k] + ) + mask = mask_counts[k] >= mask_threshold[1] + mask &= mask_counts[k] < mask_threshold[0] - if True: # limit y-axis to 1e1 - ax.set_ylim([None, 1e1]) + l = ax.plot(snr_factors[mask], time_accuracies[k][mask], **kwargs) + template_dt_colors[k] = l[0].get_color() + + # indicate threshold + if True: + ax.axhline(template_dt/np.sqrt(12), ls='--', alpha=0.7, color=template_dt_colors[k], label=f'Template dt:{template_dt:0.1e}ns') # Set horizontal line at 1 ns if True: ax.axhline(1, ls='--', alpha=0.8, color='g') - ax.grid() - ax.axhline(template_dt/np.sqrt(12), ls='--', alpha=0.7, color='b') + + ax.legend() + + if True: # limit y-axis to 1e1 + ax.set_ylim([None, 1e1]) fig.tight_layout() fig.savefig(f"figures/11_time_res_vs_snr_tdt{template_dt:0.1e}.pdf")