diff --git a/simulations/airshower_beacon_simulation/ca_period_from_shower.py b/simulations/airshower_beacon_simulation/ca_period_from_shower.py index 30435c5..1acaaac 100755 --- a/simulations/airshower_beacon_simulation/ca_period_from_shower.py +++ b/simulations/airshower_beacon_simulation/ca_period_from_shower.py @@ -189,6 +189,7 @@ if __name__ == "__main__": fname_dir = args.data_dir antennas_fname = path.join(fname_dir, beacon.antennas_fname) time_diffs_fname = 'time_diffs.hdf5' if not True else antennas_fname + tx_fname = path.join(fname_dir, beacon.tx_fname) ## This is a file indicating whether the k-finding algorithm was ## stopped early. This happens when the ks do not change between @@ -203,6 +204,7 @@ if __name__ == "__main__": # Read in antennas from file _, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) + _, __, txdata = beacon.read_tx_file(tx_fname) # Read original REvent ev = REvent(args.input_fname) # .. patch in our antennas @@ -216,53 +218,81 @@ if __name__ == "__main__": freq_name = next(iter(freq_names)) f_beacon = ev.antennas[0].beacon_info[freq_name]['freq'] - # Prepare polarisation and passbands - rit.set_pol_and_bp(ev, low=low_bp, high=high_bp) - ## ## Manipulate time and traces of each antenna ## ### Remove time due to true phase ### and optionally remove the beacon + + ### Note: there is no use in changing *_AxB variables here (except for plotting), + ### they're recomputed by the upcoming rit.set_pol_and_bp call. measured_repair_offsets = beacon.read_antenna_clock_repair_offsets(ev.antennas, mode='phases', freq_name=freq_name) for i, ant in enumerate(ev.antennas): - ev.antennas[i].orig_t = ev.antennas[i].t_AxB + ev.antennas[i].orig_t = ev.antennas[i].t + ev.antennas[i].t += measured_repair_offsets[i] + # t_AxB will be set by the rit.set_pol_and_bp function ev.antennas[i].t_AxB += measured_repair_offsets[i] if apply_signal_window_from_max: N_pre, N_post = 250, 250 # TODO: make this configurable - max_idx = np.argmax(ant.E_AxB) + # Get max idx from all the traces + # and select the strongest + max_idx = [] + maxs = [] + for trace in [ant.Ex, ant.Ey, ant.Ez]: + idx = np.argmax(np.abs(trace)) + max_idx.append(idx) + maxs.append( np.abs(trace[idx]) ) + + idx = np.argmax(maxs) + max_idx = max_idx[idx] + + # Create window around max_idx low_idx = max(0, max_idx-N_pre) high_idx = min(len(ant.t), max_idx+N_post) ev.antennas[i].orig_t = ant.orig_t[low_idx:high_idx] ev.antennas[i].t = ant.t[low_idx:high_idx] - ev.antennas[i].t_AxB = ant.t_AxB[low_idx:high_idx] ev.antennas[i].Ex = ant.Ex[low_idx:high_idx] ev.antennas[i].Ey = ant.Ey[low_idx:high_idx] ev.antennas[i].Ez = ant.Ez[low_idx:high_idx] + + ev.antennas[i].t_AxB = ant.t_AxB[low_idx:high_idx] ev.antennas[i].E_AxB = ant.E_AxB[low_idx:high_idx] - + # .. and remove the beacon from the traces + # Note: ant.E_AxB is recalculated by rit.set_pol_and_bp if remove_beacon_from_trace: clock_phase = measured_repair_offsets[i]*2*np.pi*f_beacon beacon_phase = ant.beacon_info[freq_name]['beacon_phase'] f = ant.beacon_info[freq_name]['freq'] ampl = ant.beacon_info[freq_name]['amplitude'] - calc_beacon = lib.sine_beacon(f, ev.antennas[i].t_AxB, amplitude=ampl, phase=beacon_phase-clock_phase) + calc_beacon = lib.sine_beacon(f, ev.antennas[i].t, amplitude=ampl, phase=beacon_phase-clock_phase) + tx_amps = txdata['amplitudes'] + tx_amps_sum = np.sum(tx_amps) + # Split up contribution to the various polarisations + for j, amp in enumerate(tx_amps): + if j == 0: + ev.antennas[i].Ex -= amp*(1/tx_amps_sum)*calc_beacon + elif j == 1: + ev.antennas[i].Ey -= amp*(1/tx_amps_sum)*calc_beacon + elif j == 2: + ev.antennas[i].Ez -= amp*(1/tx_amps_sum)*calc_beacon + + # Subtract the beacon from E_AxB ev.antennas[i].E_AxB -= calc_beacon # Make a figure of the manipulated traces - if i == 2: + if i == 72: orig_beacon_amplifier = ampl/max(ant.beacon) fig, ax = plt.subplots(figsize=figsize) - ax.set_title(f"Signal and Beacon traces Antenna {i}") + ax.set_title(f"Signal and Beacon traces Antenna {ant.name}") ax.set_xlabel("Time [ns]") ax.set_ylabel("Amplitude [$\\mu V/m$]") @@ -280,25 +310,27 @@ if __name__ == "__main__": old_xlim = ax.get_xlim() if True: # zoomed on part without peak of this trace - wx, x = 100, 0#ant.t_AxB[np.argmax(ant.E_AxB)] - ax.set_xlim(x-wx, x+wx) + wx, x = 200, min(ant.t_AxB) + ax.set_xlim(x-5, x+wx) - fig.savefig(path.join(fig_dir, __file__+f'.traces.A{i}.zoomed.beacon.pdf')) + fig.savefig(path.join(fig_dir, __file__+f'.traces.A{ant.name}.zoomed.beacon.pdf')) if True: # zoomed on peak of this trace idx = np.argmax(ev.antennas[i].E_AxB) x = ev.antennas[i].t_AxB[idx] - wx = 100 - ax.set_xlim(x-wx, x+wx) - fig.savefig(path.join(fig_dir, __file__+f".traces.A{i}.zoomed.peak.pdf")) + wx = 300 + ax.set_xlim(x-wx//2, x+wx//2) + fig.savefig(path.join(fig_dir, __file__+f".traces.A{ant.name}.zoomed.peak.pdf")) ax.set_xlim(*old_xlim) - fig.savefig(path.join(fig_dir, __file__+f'.traces.A{i}.pdf')) + fig.savefig(path.join(fig_dir, __file__+f'.traces.A{ant.name}.pdf')) if show_plots: plt.show() + # Prepare polarisation and passbands + rit.set_pol_and_bp(ev, low=low_bp, high=high_bp) # determine allowable ks per location dt = ev.antennas[0].t_AxB[1] - ev.antennas[0].t_AxB[0] diff --git a/simulations/airshower_beacon_simulation/da_reconstruction.py b/simulations/airshower_beacon_simulation/da_reconstruction.py index d85d764..f54ec7c 100755 --- a/simulations/airshower_beacon_simulation/da_reconstruction.py +++ b/simulations/airshower_beacon_simulation/da_reconstruction.py @@ -83,6 +83,36 @@ if __name__ == "__main__": ev.antennas[i].t += measured_repair_offsets[i] ev.antennas[i].t_AxB += measured_repair_offsets[i] + if apply_signal_window_from_max: + N_pre, N_post = 250, 250 # TODO: make this configurable + + # Get max idx from all the traces + # and select the strongest + max_idx = [] + maxs = [] + + for trace in [ant.Ex, ant.Ey, ant.Ez]: + idx = np.argmax(np.abs(trace)) + max_idx.append(idx) + maxs.append( np.abs(trace[idx]) ) + + idx = np.argmax(maxs) + max_idx = max_idx[idx] + + # Create window around max_idx + low_idx = max(0, max_idx-N_pre) + high_idx = min(len(ant.t), max_idx+N_post) + + ev.antennas[i].orig_t = ant.orig_t[low_idx:high_idx] + ev.antennas[i].t = ant.t[low_idx:high_idx] + + ev.antennas[i].Ex = ant.Ex[low_idx:high_idx] + ev.antennas[i].Ey = ant.Ey[low_idx:high_idx] + ev.antennas[i].Ez = ant.Ez[low_idx:high_idx] + + ev.antennas[i].t_AxB = ant.t_AxB[low_idx:high_idx] + ev.antennas[i].E_AxB = ant.E_AxB[low_idx:high_idx] + # .. and remove the beacon from the traces # Note: ant.E_AxB is recalculated by rit.set_pol_and_bp if remove_beacon_from_traces: @@ -111,7 +141,7 @@ if __name__ == "__main__": ## ## Make a figure of the manipulated traces ## - if i == 2: + if i == 72: orig_beacon_amplifier = ampl_AxB/max(ant.beacon) for k in range(2): @@ -127,7 +157,7 @@ if __name__ == "__main__": fname_extra = ".Ex" fig, ax = plt.subplots(figsize=figsize) - ax.set_title(f"Signal and Beacon traces Antenna {i}") + ax.set_title(f"Signal and Beacon traces Antenna {ant.name}") ax.set_xlabel("Time [ns]") ax.set_ylabel("Amplitude [$\\mu V/m$]") @@ -151,14 +181,14 @@ if __name__ == "__main__": wx, x = 100, 100 ax.set_xlim(x-wx, x+wx) - fig.savefig(path.join(fig_dir, __file__+f'.traces.A{i}.zoomed.beacon{fname_extra}.pdf')) + fig.savefig(path.join(fig_dir, __file__+f'.traces.A{ant.name}.zoomed.beacon{fname_extra}.pdf')) if True: # zoomed on peak of this trace idx = np.argmax(ev.antennas[i].E_AxB) x = ev.antennas[i].t_AxB[idx] wx = 100 ax.set_xlim(x-wx, x+wx) - fig.savefig(path.join(fig_dir, __file__+f".traces.A{i}.zoomed.peak{fname_extra}.pdf")) + fig.savefig(path.join(fig_dir, __file__+f".traces.A{ant.name}.zoomed.peak{fname_extra}.pdf")) ax.set_xlim(*old_xlim) diff --git a/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py b/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py index 210e158..d1e460d 100755 --- a/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py +++ b/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py @@ -135,7 +135,7 @@ if __name__ == "__main__": elif j == 2: ev.antennas[i].Ez -= amp*(1/tx_amps_sum)*calc_beacon - # + # Subtract the beacon from E_AxB ev.antennas[i].E_AxB -= calc_beacon # Slice the traces to a small part around the peak @@ -143,7 +143,19 @@ if __name__ == "__main__": N_pre, N_post = 250, 250 # TODO: make this configurable for i, ant in enumerate(ev.antennas): - max_idx = np.argmax(ant.E_AxB) + # Get max idx from all the traces + # and select the strongest + max_idx = [] + maxs = [] + + for trace in [ant.Ex, ant.Ey, ant.Ez]: + idx = np.argmax(np.abs(trace)) + max_idx.append(idx) + maxs.append( np.abs(trace[idx]) ) + + idx = np.argmax(maxs) + max_idx = max_idx[idx] + low_idx = max(0, max_idx-N_pre) high_idx = min(len(ant.t), max_idx+N_post) @@ -156,13 +168,13 @@ if __name__ == "__main__": ev.antennas[i].Ez = ant.Ez[low_idx:high_idx] ev.antennas[i].E_AxB = ant.E_AxB[low_idx:high_idx] + ## Apply polarisation and bandpass filter + rit.set_pol_and_bp(ev) + # backup antenna times backup_antenna_t = [ ant.t for ant in ev.antennas ] backup_antenna_t_AxB = [ ant.t_AxB for ant in ev.antennas ] - ## Apply polarisation and bandpass filter - rit.set_pol_and_bp(ev) - with joblib.parallel_backend("loky"): for case in wanted_cases: print(f"Starting {case} figure")