diff --git a/simulations/11_pulsed_timing.py b/simulations/11_pulsed_timing.py index bdabee5..b8fa371 100755 --- a/simulations/11_pulsed_timing.py +++ b/simulations/11_pulsed_timing.py @@ -160,6 +160,18 @@ def write_time_residuals_cache(cache_fname, time_residuals, template_dt, antenna ds = pgroup2.create_dataset(ds_name, (len(time_residuals)), dtype='f', data=time_residuals, maxshape=(None)) +def create_template(dt=1, timelength=1, bp_freq=(0, np.inf), name=None, normalise=False): + template = Waveform(None, dt=dt, name=name) + _deltapeak = util.deltapeak(timelength=timelength, samplerate=1/dt, offset=0) + template.signal = antenna_bp(_deltapeak[0], *bp_freq, dt) + template.peak_sample = _deltapeak[1] + template.peak_time = template.dt * template.peak_sample + + if normalise: + template.signal /= max(template.signal) + + return template, _deltapeak + if __name__ == "__main__": import os import matplotlib @@ -186,36 +198,26 @@ if __name__ == "__main__": antenna_timelength = 1024 # ns cut_wrong_peak_matches = True + normalise_noise = False + + h5_cache_fname = f'11_pulsed_timing.hdf5' # - # Create the template - # - template = Waveform(None, dt=template_dt, name='Template') - _deltapeak = util.deltapeak(timelength=template_length, samplerate=1/template.dt, offset=0) - template.signal = antenna_bp(_deltapeak[0], *bp_freq, template_dt) - template.peak_sample = _deltapeak[1] - template.peak_time = template.dt * template.peak_sample - # Interpolation Template # to create an 'analog' sampled antenna - interp_template = Waveform(None, dt=interp_template_dt, name='Interpolation Template') - _interp_deltapeak = util.deltapeak(timelength=template_length, samplerate=1/interp_template.dt, offset=0) - interp_template.signal = antenna_bp(_interp_deltapeak[0], *bp_freq, interp_template.dt) - interp_template.peak_sample = _interp_deltapeak[1] - interp_template.peak_time = interp_template.dt * interp_template.peak_sample - interp_template.signal *= max(template.signal)/max(interp_template.signal) - interp_template.interpolate = interpolate.interp1d(interp_template.t, interp_template.signal, assume_sorted=True, fill_value=0, bounds_error=False) + interp_template, _deltapeak = create_template(dt=interp_template_dt, timelength=template_length, bp_freq=bp_freq, name='Interpolation Template', normalise=True) - if True: # show template + interp_template.interpolate = interpolate.interp1d(interp_template.t, interp_template.signal, assume_sorted=True, fill_value=0, bounds_error=False, copy=False) + + if True: # show interpolation template fig, ax = plt.subplots() ax.set_title("Deltapeak and Bandpassed Template") ax.set_xlabel("Time [ns]") ax.set_ylabel("Amplitude") - ax.plot(template.t, max(template.signal)*_deltapeak[0], label='Impulse Template') - ax.plot(template.t, template.signal, label='Filtered Template') - ax.plot(interp_template.t, interp_template.signal, label='Filtered Interpolation Template', ls='dashed') + ax.plot(interp_template.t, max(interp_template.signal)*_deltapeak[0], label='Impulse Template') + ax.plot(interp_template.t, interp_template.signal, label='Filtered Template') ax.legend() - fig.savefig('figures/11_template_deltapeak.pdf') + fig.savefig('figures/11_interpolation_deltapeak.pdf') if True: # show filtering equivalence samplerates _deltapeak = util.deltapeak(timelength=template_length, samplerate=1/antenna_dt, offset=0) @@ -226,16 +228,20 @@ if __name__ == "__main__": ax.plot(_time, _bandpassed, label='Filtered Antenna') ax.legend() - fig.savefig('figures/11_template_deltapeak+antenna.pdf') + fig.savefig('figures/11_interpolation_deltapeak+antenna.pdf') 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 # - 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)): @@ -256,7 +262,7 @@ if __name__ == "__main__": ## place the deltapeak signal at a random location antenna = Waveform(None, dt=antenna_dt, name='Signal') - if False: # Create antenna trace without template + if False: # Create antenna trace without interpolation template antenna_true_signal, antenna_peak_sample = util.deltapeak(timelength=antenna_timelength, samplerate=1/antenna.dt, offset=[0.2, 0.8], rng=rng) antenna.peak_sample = antenna_peak_sample @@ -265,7 +271,7 @@ if __name__ == "__main__": print(f"Antenna Peak Time: {antenna.peak_time}") print(f"Antenna Peak Sample: {antenna.peak_sample}") - else: # Sample the template at some offset + else: # Sample the interpolation template at some offset antenna.peak_time = antenna_timelength * ((0.8 - 0.2) *rng.random(1) + 0.2) sampling_offset = rng.random(1)*antenna.dt @@ -284,7 +290,7 @@ if __name__ == "__main__": ## Add noise noise_amplitude = max(template.signal) * 1/snr_sigma_factor - noise_realisation = noise_amplitude * white_noise_realisation(len(antenna.signal), normalise=False) + noise_realisation = noise_amplitude * white_noise_realisation(len(antenna.signal), normalise=normalise_noise) filtered_noise = antenna_bp(noise_realisation, *bp_freq, antenna.dt) antenna.signal += filtered_noise @@ -569,7 +575,7 @@ if __name__ == "__main__": l = ax.plot(snr_factors[mask], time_accuracies[mask], **kwargs) - if True: # limit y-axis to 1e0 + if True: # limit y-axis to 1e1 ax.set_ylim([None, 1e1])