mirror of
https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git
synced 2024-11-13 18:13:31 +01:00
102 lines
2.9 KiB
Python
Executable file
102 lines
2.9 KiB
Python
Executable file
#!/usr/bin/env python3
|
|
# vim: fdm=indent ts=4
|
|
|
|
import h5py
|
|
import matplotlib.pyplot as plt
|
|
from mpl_toolkits.mplot3d import Axes3D
|
|
import numpy as np
|
|
from copy import deepcopy
|
|
|
|
import aa_generate_beacon as beacon
|
|
import lib
|
|
from lib import rit
|
|
from earsim import REvent, Antenna
|
|
from atmocal import AtmoCal
|
|
import pickle
|
|
|
|
if __name__ == "__main__":
|
|
from os import path
|
|
import sys
|
|
|
|
fname = "ZH_airshower/mysim.sry"
|
|
|
|
####
|
|
fname_dir = path.dirname(fname)
|
|
if True: # use the original files
|
|
ev = REvent(fname)
|
|
|
|
# Full reconstruction?
|
|
if False:
|
|
print("Full reconstruction")
|
|
with open('ritdir/res.pkl', 'wb') as fp:
|
|
res = rit.reconstruction(ev, outfile='ritdir/fig.pdf', slice_outdir='ritdir/')
|
|
|
|
pickle.dump(res, fp)
|
|
sys.exit(0)
|
|
|
|
# Else do single slice for original timing
|
|
# and randomised timing
|
|
print("Shower plane reconstruction + randomised timing")
|
|
|
|
atm = AtmoCal()
|
|
X = 750
|
|
max_clock_offset = 10# ns
|
|
rit.set_pol_and_bp(ev)
|
|
dXref = atm.distance_to_slant_depth(np.deg2rad(ev.zenith),X,0)
|
|
scale2d = dXref*np.tan(np.deg2rad(2.))
|
|
|
|
ants_copy = deepcopy(ev.antennas)
|
|
# randomise timing for antenna copies
|
|
seed = 12345
|
|
rng = np.random.default_rng(seed)
|
|
|
|
clock_offsets = max_clock_offset * (2*rng.uniform(size=len(ants_copy)) - 1)
|
|
for i, ant in enumerate(ants_copy):
|
|
ant.t -= clock_offsets[i]
|
|
|
|
fig, axs = plt.subplots(1,2, sharex=True, sharey=True)
|
|
fig.suptitle(f"Reconstructed Shower plane slice at X={X}")
|
|
for ax in axs:
|
|
ax.set_xlabel("-v x B (m)")
|
|
ax.set_ylabel(" vxvxB (m)")
|
|
axs[0].set_title("Simulated timing")
|
|
axs[1].set_title(f"Randomised timing [0, {max_clock_offset}]")
|
|
|
|
try:
|
|
for i, ax in enumerate(axs):
|
|
# modify timing
|
|
if i == 1:
|
|
ev.antennas = ants_copy
|
|
|
|
xx,yy,p,km= rit.shower_plane_slice(ev,X,21,21,scale2d,scale2d)
|
|
|
|
if i == 0:
|
|
vmax, vmin = max(p), min(p)
|
|
|
|
sc = ax.scatter(xx,yy,c=p, vmax=vmax, vmin=vmin)
|
|
|
|
# Plot centers
|
|
im = np.argmax(p)
|
|
|
|
ax.plot(0, 0,'r+',ms=30)
|
|
ax.plot(xx[im],yy[im],'bx',ms=30)
|
|
|
|
fig.colorbar(sc, ax=axs[0])
|
|
# Always make sure to show the plot
|
|
except:
|
|
plt.show()
|
|
|
|
if True:
|
|
fig.savefig(f"shower_plane_slice_X{X}_T{max_clock_offset}.pdf")
|
|
|
|
|
|
else:
|
|
raise NotImplementedError
|
|
antennas_fname = path.join(fname_dir, beacon.antennas_fname)
|
|
|
|
if not path.isfile(antennas_fname):
|
|
print("Antenna file cannot be found, did you try generating a beacon?")
|
|
sys.exit(1)
|
|
|
|
|
|
plt.show()
|