m-thesis-introduction/simulations/airshower_beacon_simulation/da_reconstruction.py

#!/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
from mpl_toolkits.mplot3d import Axes3D # required for projection='3d' on old matplotliblib versions
import numpy as np
from os import path
import pickle
import joblib

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()

    from scriptlib import MyArgumentParser
    parser = MyArgumentParser()
    args = parser.parse_args()

    fname = "ZH_airshower/mysim.sry"

    fig_dir = args.fig_dir
    fig_subdir = path.join(fig_dir, 'reconstruction')
    show_plots = args.show_plots

    ####
    fname_dir = path.dirname(fname)
    antennas_fname = path.join(fname_dir, beacon.antennas_fname)
    pickle_fname = path.join(fname_dir, 'res.pkl')

    # create fig_dir
    if fig_dir:
        os.makedirs(fig_dir, exist_ok=True)
        if fig_subdir:
            os.makedirs(fig_subdir, exist_ok=True)

    # 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

    # For now only implement using one freq_name
    freq_names = antennas[0].beacon_info.keys()
    if len(freq_names) > 1:
        raise NotImplementedError

    freq_name = next(iter(freq_names))
    f_beacon = ev.antennas[0].beacon_info[freq_name]['freq']

    # Repair clock offsets with the measured offsets
    measured_repair_offsets = beacon.read_antenna_clock_repair_offsets(ev.antennas, mode='all', freq_name=freq_name)
    for i, ant in enumerate(ev.antennas):
        # t_AxB will be set by the rit.set_pol_and_bp function
        ev.antennas[i].t += measured_repair_offsets[i]

    N_X, Xlow, Xhigh = 23, 100, 1200
    with joblib.parallel_backend("loky"):
        res = rit.reconstruction(ev, outfile=fig_subdir+'/fig.pdf', slice_outdir=fig_subdir+'/', Xlow=Xlow, N_X=N_X, Xhigh=Xhigh)

    ## Save a pickle
    with open(pickle_fname, 'wb') as fp:
        pickle.dump(res,fp)

    if show_plots:
        plt.show()