ZH: shower slice figures for best k set per iteration

Extends the reconstruction from 49d4779
This commit is contained in:
Eric Teunis de Boone 2023-01-11 02:20:06 +01:00
parent d5d1686a6b
commit 32f7445fc4

View file

@ -339,19 +339,23 @@ if __name__ == "__main__":
if True: #plot maximum at test locations if True: #plot maximum at test locations
fig, axs = plt.subplots() fig, axs = plt.subplots()
axs.set_title(f"Optimizing signal strength, Grid Run {r}") axs.set_title(f"Optimizing signal strength varying k per antenna,\n Grid Run {r}")
axs.set_ylabel("vxvxB [km]") axs.set_ylabel("vxvxB [km]")
axs.set_xlabel(" vxB [km]") axs.set_xlabel(" vxB [km]")
axs.set_aspect('equal', 'datalim') axs.set_aspect('equal', 'datalim')
sc = axs.scatter(xx/1e3, yy/1e3, c=maxima_per_loc, cmap='Spectral_r', alpha=0.6) sc = axs.scatter(xx/1e3, yy/1e3, c=maxima_per_loc, cmap='Spectral_r', alpha=0.6)
fig.colorbar(sc, ax=axs) fig.colorbar(sc, ax=axs, label='Max Amplitude [$\\mu V/m$]')
# indicate maximum value
idx = np.argmax(maxima_per_loc)
axs.plot(xx[idx]/1e3, yy[idx]/1e3, 'bx', ms=30)
if fig_dir: if fig_dir:
old_xlims = axs.get_xlim() old_xlims = axs.get_xlim()
old_ylims = axs.get_ylim() old_ylims = axs.get_ylim()
fig.tight_layout() fig.tight_layout()
fig.savefig(path.join(fig_dir, __file__+f'.maxima.run{r}.pdf')) fig.savefig(path.join(fig_dir, __file__+f'.maxima.run{r}.pdf'))
if True: if False:
axs.plot(tx.x/1e3, tx.y/1e3, marker='X', color='k') axs.plot(tx.x/1e3, tx.y/1e3, marker='X', color='k')
fig.tight_layout() fig.tight_layout()
fig.savefig(path.join(fig_dir, __file__+f'.maxima.run{r}.with_tx.pdf')) fig.savefig(path.join(fig_dir, __file__+f'.maxima.run{r}.with_tx.pdf'))
@ -372,19 +376,51 @@ if __name__ == "__main__":
## Do a small reconstruction of the shower for best ks ## Do a small reconstruction of the shower for best ks
if True: if True:
print("Reconstructing for best k") print("Reconstructing for best k")
_, __, p, ___ = rit.shower_plane_slice(ev, X=Xref, Nx=len(x), Ny=len(y), wx=scale2d, wy=scale2d, xoff=xoff, yoff=yoff, zgr=0)
for j in range(2):
power_reconstruction = j==1
if power_reconstruction: # Do power reconstruction
# backup antenna times
backup_times = [ ant.t_AxB for ant in ev.antennas ]
# incorporate ks into timing
for i, ant in enumerate(ev.antennas):
ev.antennas[i].t_AxB = ant.t_AxB - best_k[i] * 1/f_beacon
xx, yy, p, ___ = rit.shower_plane_slice(ev, X=Xref, Nx=len(x), Ny=len(y), wx=x[-1]-x[0], wy=y[-1]-y[0], xoff=xoff, yoff=yoff, zgr=0)
# repair antenna times
for i, backup_t_AxB in enumerate(backup_times):
ev.antennas[i].t_AxB = backup_t_AxB
else: # get maximum amplitude at each location
maxima = np.empty( len(locs) )
for i, loc in enumerate(locs):
test_loc = loc[0]* ev.uAxB + loc[1]*ev.uAxAxB + dXref *ev.uA
P, t_, a_, a_sum, t_sum = rit.pow_and_time(test_loc, ev, dt=dt)
maxima[i] = np.max(a_sum)
fig, axs = plt.subplots() fig, axs = plt.subplots()
axs.set_title(f"Shower slice for best k, Grid Run {r}") axs.set_title(f"Shower slice for best k, Grid Run {r}")
axs.set_ylabel("vxvxB [km]") axs.set_ylabel("vxvxB [km]")
axs.set_xlabel(" vxB [km]") axs.set_xlabel(" vxB [km]")
axs.set_aspect('equal', 'datalim') axs.set_aspect('equal', 'datalim')
if power_reconstruction:
sc = axs.scatter(xx/1e3, yy/1e3, c=p, cmap='Spectral_r', alpha=0.6) sc = axs.scatter(xx/1e3, yy/1e3, c=p, cmap='Spectral_r', alpha=0.6)
fig.colorbar(sc, ax=axs) fig.colorbar(sc, ax=axs, label='Power')
else:
sc = axs.scatter(xx/1e3, yy/1e3, c=maxima, cmap='Spectral_r', alpha=0.6)
fig.colorbar(sc, ax=axs, label='Max Amplitude [$\\mu V/m$]')
if fig_dir: if fig_dir:
if power_reconstruction:
fname_extra = "power"
else:
fname_extra = "max_amp"
fig.tight_layout() fig.tight_layout()
fig.savefig(path.join(fig_dir, __file__+f'.reconstruction.run{r}.pdf')) fig.savefig(path.join(fig_dir, __file__+f'.reconstruction.run{r}.{fname_extra}.pdf'))
# Abort if no improvement # Abort if no improvement
if ( r!= 0 and (old_ks_per_loc == ks_per_loc[best_idx]).all() ): if ( r!= 0 and (old_ks_per_loc == ks_per_loc[best_idx]).all() ):