ZH: reconstruction script update

simulations/airshower_beacon_simulation/do_reconstruction.py

@@ -3,7 +3,7 @@
 
 import h5py
 import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
+from mpl_toolkits.mplot3d import Axes3D # required for projection='3d' on old matplotliblib versions
 import numpy as np
 from copy import deepcopy
 
@@ -19,17 +19,35 @@ if __name__ == "__main__":
     import sys
 
     fname = "ZH_airshower/mysim.sry"
+    dirname = '/scratch/etdeboone/reconstructions/'+  (sys.argv[2]+'/' if 2 in sys.argv else '')
+    print("Dirname:", dirname)
 
     ####
     fname_dir = path.dirname(fname)
     if True: # use the original files
         ev = REvent(fname)
 
+        max_clock_offset = int(sys.argv[1]) if len(sys.argv) > 1 else 5# ns
+        # randomise timing for antenna copies
+        seed = 12345
+        rng = np.random.default_rng(seed)
+
+        dirname += f'rand{max_clock_offset}ns'
+
         # 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/')
+
+            if max_clock_offset:
+                print(f"Randomising timing: {max_clock_offset}")
+                clock_offsets = max_clock_offset * (2*rng.uniform(size=len(ev.antennas)) - 1)
+                for i, ant in enumerate(ev.antennas):
+                    ant.t -= clock_offsets[i]
+            else:
+                print(f"Not randomising timing")
+
+            with open(dirname + '/res.pkl', 'wb') as fp:
+                res = rit.reconstruction(ev, outfile=dirname+'/fig.pdf', slice_outdir=dirname+'/', Xlow=100, N_X=15, Xhigh=800)
 
                 pickle.dump(res, fp)
             sys.exit(0)
@@ -40,16 +58,11 @@ if __name__ == "__main__":
 
         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]
@@ -87,7 +100,7 @@ if __name__ == "__main__":
             plt.show()
 
         if True:
-            fig.savefig(f"shower_plane_slice_X{X}_T{max_clock_offset}.pdf")
+            fig.savefig(dirname + f"/shower_plane_slice_X{X}_T{max_clock_offset}.pdf")
             
 
     else:

simulations/airshower_beacon_simulation/lib/rit.py

@@ -318,7 +318,7 @@ def fill_stations_propeties(e,res):
     res.station_ids.append(ids)
     #res.has_pulse.append(has_pulse)
 
-def reconstruction(e,outfile='', slice_outdir=None):
+def reconstruction(e,outfile='', slice_outdir=None, Xlow=300, Xhigh=1000, N_X=15):
     res = RITResult()
     res.isMC.append(True)
     res.zenith_ini.append(e.zenith)
@@ -329,7 +329,7 @@ def reconstruction(e,outfile='', slice_outdir=None):
 
     #only use signal that have a signal in data
     fill_stations_propeties(e,res)
-    Xs,axis_points,max_vals = get_axis_points(e,slice_outdir is not None, slice_outdir, 200, 1000, 5)
+    Xs,axis_points,max_vals = get_axis_points(e,slice_outdir is not None, slice_outdir, Xlow, Xhigh, N_X)
     zen,azi,core = fit_track(e,axis_points,max_vals,1e2)
     fig = figure_3D(axis_points,max_vals,zen,azi,core,res)
     fig.savefig(outfile)