From ffb1fa8e42578aa206e062d0650aef596f3ae36b Mon Sep 17 00:00:00 2001
From: Eric Teunis de Boone <ericteunis@deboone.nl>
Date: Mon, 9 Jan 2023 16:40:32 +0100
Subject: [PATCH] ZH: init cc full reconstruction script after fixing clock
 offsets

---
 .../cc_power_along_axis_reconstruction.py     | 63 +++++++++++++++++++
 1 file changed, 63 insertions(+)
 create mode 100755 simulations/airshower_beacon_simulation/cc_power_along_axis_reconstruction.py

diff --git a/simulations/airshower_beacon_simulation/cc_power_along_axis_reconstruction.py b/simulations/airshower_beacon_simulation/cc_power_along_axis_reconstruction.py
new file mode 100755
index 0000000..0e4399b
--- /dev/null
+++ b/simulations/airshower_beacon_simulation/cc_power_along_axis_reconstruction.py
@@ -0,0 +1,63 @@
+#!/usr/bin/env python3
+# vim: fdm=indent ts=4
+
+"""
+Do a reconstruction of airshower after correcting for the
+clock offsets.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+from os import path
+
+from earsim import REvent
+from atmocal import AtmoCal
+import aa_generate_beacon as beacon
+import lib
+from lib import rit
+
+
+if __name__ == "__main__":
+    import sys
+    import os
+    import matplotlib
+    if os.name == 'posix' and "DISPLAY" not in os.environ:
+        matplotlib.use('Agg')
+
+    atm = AtmoCal()
+
+    fname = "ZH_airshower/mysim.sry"
+
+    fig_dir = "./figures/"
+    fig_subdir = path.join(fig_dir, 'reconstruction')
+    show_plots = False
+
+    ####
+    fname_dir = path.dirname(fname)
+    antennas_fname = path.join(fname_dir, beacon.antennas_fname)
+    time_diffs_fname = 'time_diffs.hdf5' if not True else antennas_fname
+
+    # Read in antennas from file
+    _, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
+    # Read original REvent
+    ev = REvent(fname)
+    # .. patch in our antennas
+    ev.antennas = antennas
+
+    # Repair clock offsets with the measured offsets
+    for i, ant in enumerate(ev.antennas):
+        clock_phase_time = ant.beacon_info[freq_name]['sigma_phase_mean']/(2*np.pi*f_beacon)
+        best_k_time = ant.beacon_info[freq_name]['best_k_time']
+
+        total_clock_time = best_k_time + clock_phase_time
+        ev.antennas[i].t += total_clock_time
+
+    ## Start a pickle
+    with open(pickle_fname, 'wb') as fp:
+        N_X, Xlow, Xhigh = 23, 100, 1200
+        res = rit.reconstruction(ev, outfile=fig_subdir+'/fig.pdf', slice_outdir=fig_subdir, Xlow=Xlow, N_X=N_X, Xhigh=Xhigh)
+
+        pickle.dump(res,fp)
+
+    if show_plots():
+        plt.show()