From b71676574531563103702ce30bec1fe407bfe88b Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Tue, 29 Nov 2022 17:04:26 +0100 Subject: [PATCH] ZH: rename period shower script --- .../bb_beacon_true_phases.py | 2 +- ...ower_plane.py => bc_period_from_shower.py} | 81 +++++++++++++------ 2 files changed, 56 insertions(+), 27 deletions(-) rename simulations/airshower_beacon_simulation/{bc_period_shower_plane.py => bc_period_from_shower.py} (75%) diff --git a/simulations/airshower_beacon_simulation/bb_beacon_true_phases.py b/simulations/airshower_beacon_simulation/bb_beacon_true_phases.py index 1127158..f826a20 100755 --- a/simulations/airshower_beacon_simulation/bb_beacon_true_phases.py +++ b/simulations/airshower_beacon_simulation/bb_beacon_true_phases.py @@ -28,7 +28,7 @@ def antenna_true_phases(tx, antennas, freq_name, c_light=3e8): geom_time = lib.geometry_time(tx, ant, c_light=c_light) geom_phase = geom_time * 2*np.pi*f_beacon - true_phases[i] = lib.phase_mod( lib.phase_mod(geom_phase) - lib.phase_mod(measured_phase)) + true_phases[i] = lib.phase_mod(lib.phase_mod(measured_phase) - lib.phase_mod(geom_phase) ) return true_phases diff --git a/simulations/airshower_beacon_simulation/bc_period_shower_plane.py b/simulations/airshower_beacon_simulation/bc_period_from_shower.py similarity index 75% rename from simulations/airshower_beacon_simulation/bc_period_shower_plane.py rename to simulations/airshower_beacon_simulation/bc_period_from_shower.py index 4e86d3d..8d2a434 100755 --- a/simulations/airshower_beacon_simulation/bc_period_shower_plane.py +++ b/simulations/airshower_beacon_simulation/bc_period_from_shower.py @@ -3,7 +3,7 @@ """ Find the best integer multiple to shift -antennas to be able to resolve +antennas to be able to resolve """ import h5py @@ -18,7 +18,7 @@ import aa_generate_beacon as beacon import lib from lib import rit -def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_sample_shifts, dt=1): +def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_sample_shifts, dt=None, sample_shift_first_trace=0): """ Find the best sample_shift for each antenna by summing the antenna traces and seeing how to get the best alignment. @@ -28,6 +28,9 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp t_min = 1e9 t_max = -1e9 + if dt is None: + dt = antennas[0].t_AxB[1] - antennas[0].t_AxB[0] + # propagate to test location for i, ant in enumerate(antennas): aloc = [ant.x, ant.y, ant.z] @@ -52,13 +55,12 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp best_sample_shifts = np.zeros( (len(antennas)) ,dtype=int) for i, (t_r, E_) in enumerate(zip(t_, a_)): - t_min = t_r[0] f = interp1d(t_r, E_, assume_sorted=True, bounds_error=False, fill_value=0) a_int = f(t_sum) if i == 0: a_sum += a_int - best_sample_shifts[i] = 0 + best_sample_shifts[i] = sample_shift_first_trace continue shift_maxima = np.zeros( len(allowed_sample_shifts) ) @@ -84,10 +86,13 @@ if __name__ == "__main__": fname = "ZH_airshower/mysim.sry" allowed_ks = np.arange(-2, 3, 1) - Xref = 450 + Xref = 400 - x_coarse = np.linspace(-2e3, 2e3, 11) - y_coarse = np.linspace(-2e3, 2e3, 11) + x_coarse = np.linspace(-20e3, 20e3, 10) + y_coarse = np.linspace(-20e3, 20e3, 10) + + x_fine = np.linspace(-2e3, 2e3, 30) + y_fine = np.linspace(-2e3, 2e3, 30) #### fname_dir = path.dirname(fname) @@ -111,7 +116,8 @@ if __name__ == "__main__": # determine best ks per location dt = ev.antennas[0].t[1] - ev.antennas[0].t[0] - allowed_sample_shifts = np.ceil(allowed_ks/f_beacon /dt).astype(int) + allowed_sample_shifts = np.rint(allowed_ks/f_beacon /dt).astype(int) + print("Checking:", allowed_ks, ": shifts :", allowed_sample_shifts) # Remove time due to true phase @@ -134,17 +140,18 @@ if __name__ == "__main__": x = x_coarse y = y_coarse else: - best_idx = np.argmax(maxima_per_loc) - xoff = xx[best_idx] - yoff = yy[best_idx] - x /= 4 - y /= 4 + if r == 1: + x = x_fine + y = y_fine + else: + x /= 4 + y /= 4 - print(f"Testing grid {r} centered on {xoff}, {yoff}") + print(f"Testing grid run {r} centered on {xoff}, {yoff}") - ks_per_loc = np.zeros( (len(x)*len(y), len(ev.antennas)) ) + ks_per_loc = np.zeros( (len(x)*len(y), len(ev.antennas)) , dtype=int) maxima_per_loc = np.zeros( (len(x)*len(y))) - + ## Check each location on grid xx = [] yy = [] @@ -155,23 +162,18 @@ if __name__ == "__main__": test_loc = (x_+xoff)* ev.uAxB + (y_+yoff)*ev.uAxAxB + dXref *ev.uA xx.append(x_+xoff) yy.append(y_+yoff) - + # Find best k for each antenna shifts, maximum = find_best_sample_shifts_summing_at_location(test_loc, ev.antennas, allowed_sample_shifts, dt) - + # Translate sample shifts back into period multiple k - ks = shifts*f_beacon*dt - + ks = np.rint(shifts*f_beacon*dt) + ks_per_loc[i] = ks maxima_per_loc[i] = maximum - + xx = np.array(xx) yy = np.array(yy) - - best_idx = np.argmax(maxima_per_loc) - np.savetxt(__file__+f'.run{r}.txt') - print(ks_per_loc[best_idx]) - if True: #plot maximum at test locations fig, axs = plt.subplots() axs.set_title(f"Grid Run {r}") @@ -181,4 +183,31 @@ if __name__ == "__main__": fig.colorbar(sc, ax=axs) fig.savefig(__file__+f'.run{r}.pdf') + + ## Save ks to file + best_idx = np.argmax(maxima_per_loc) + np.savetxt(__file__+f'.bestk.run{r}.txt', ks_per_loc[best_idx] ) + print(ks_per_loc[best_idx]) + + ## Save maxima to file + np.savetxt(__file__+f'.maxima.run{r}.txt', maxima_per_loc) + + # Abort if no improvement + if ( r!= 0 and (old_ks_per_loc == ks_per_loc[best_idx]).all() ): + print("No improvement, breaking") + break + + # Prepare variables for next loop + old_ks_per_loc = ks_per_loc[best_idx] + xoff = xx[best_idx] + yoff = yy[best_idx] + + # Save best ks to hdf5 antenna file + with h5py.File(antennas_fname, 'a') as fp: + group = fp.require_group('antennas') + + for i, ant in enumerate(antennas): + h5ant = group[ant.name] + h5ant.attrs['best_k'] = old_ks_per_loc[i] + plt.show()