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

#!/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

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")
            res = rit.reconstruction(ev, outfile='ritdir/fig.pdf')
            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 = 100 # 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
        finally:
            plt.show()

            

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