ZH: rename period shower script

simulations/airshower_beacon_simulation/bb_beacon_true_phases.py

@@ -28,7 +28,7 @@ def antenna_true_phases(tx, antennas, freq_name, c_light=3e8):
         geom_time = lib.geometry_time(tx, ant, c_light=c_light)
         geom_phase = geom_time * 2*np.pi*f_beacon
 
-        true_phases[i] =  lib.phase_mod( lib.phase_mod(geom_phase) - lib.phase_mod(measured_phase))
+        true_phases[i] =  lib.phase_mod(lib.phase_mod(measured_phase) - lib.phase_mod(geom_phase) )
 
     return true_phases
 

simulations/airshower_beacon_simulation/bc_period_from_shower.py

@@ -18,7 +18,7 @@ import aa_generate_beacon as beacon
 import lib
 from lib import rit
 
-def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_sample_shifts, dt=1):
+def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_sample_shifts, dt=None, sample_shift_first_trace=0):
     """
     Find the best sample_shift for each antenna by summing the antenna traces
     and seeing how to get the best alignment.
@@ -28,6 +28,9 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp
     t_min = 1e9
     t_max = -1e9
 
+    if dt is None:
+        dt = antennas[0].t_AxB[1] - antennas[0].t_AxB[0]
+
     # propagate to test location
     for i, ant in enumerate(antennas):
         aloc = [ant.x, ant.y, ant.z]
@@ -52,13 +55,12 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp
 
     best_sample_shifts = np.zeros( (len(antennas)) ,dtype=int)
     for i, (t_r, E_) in enumerate(zip(t_, a_)):
-        t_min = t_r[0]
         f = interp1d(t_r, E_, assume_sorted=True, bounds_error=False, fill_value=0)
         a_int = f(t_sum)
 
         if i == 0:
             a_sum += a_int
-            best_sample_shifts[i] = 0
+            best_sample_shifts[i] = sample_shift_first_trace
             continue
 
         shift_maxima = np.zeros( len(allowed_sample_shifts) )
@@ -84,10 +86,13 @@ if __name__ == "__main__":
     fname = "ZH_airshower/mysim.sry"
 
     allowed_ks = np.arange(-2, 3, 1)
-    Xref = 450
+    Xref = 400
 
-    x_coarse = np.linspace(-2e3, 2e3, 11)
+    x_coarse = np.linspace(-20e3, 20e3, 10)
-    y_coarse = np.linspace(-2e3, 2e3, 11)
+    y_coarse = np.linspace(-20e3, 20e3, 10)
+
+    x_fine = np.linspace(-2e3, 2e3, 30)
+    y_fine = np.linspace(-2e3, 2e3, 30)
 
     ####
     fname_dir = path.dirname(fname)
@@ -111,7 +116,8 @@ if __name__ == "__main__":
 
     # determine best ks per location
     dt = ev.antennas[0].t[1] - ev.antennas[0].t[0]
-    allowed_sample_shifts = np.ceil(allowed_ks/f_beacon /dt).astype(int)
+    allowed_sample_shifts = np.rint(allowed_ks/f_beacon /dt).astype(int)
+    print("Checking:", allowed_ks, ": shifts :", allowed_sample_shifts)
 
 
     # Remove time due to true phase
@@ -134,15 +140,16 @@ if __name__ == "__main__":
             x = x_coarse
             y = y_coarse
         else:
-            best_idx = np.argmax(maxima_per_loc)
+            if r == 1:
-            xoff = xx[best_idx]
+                x = x_fine
-            yoff = yy[best_idx]
+                y = y_fine
-            x /= 4
+            else:
-            y /= 4
+                x /= 4
+                y /= 4
 
-        print(f"Testing grid {r} centered on {xoff}, {yoff}")
+        print(f"Testing grid run {r} centered on {xoff}, {yoff}")
 
-        ks_per_loc = np.zeros( (len(x)*len(y), len(ev.antennas)) )
+        ks_per_loc = np.zeros( (len(x)*len(y), len(ev.antennas)) , dtype=int)
         maxima_per_loc = np.zeros( (len(x)*len(y)))
 
         ## Check each location on grid
@@ -160,18 +167,13 @@ if __name__ == "__main__":
             shifts, maximum = find_best_sample_shifts_summing_at_location(test_loc, ev.antennas, allowed_sample_shifts, dt)
 
             # Translate sample shifts back into period multiple k
-            ks = shifts*f_beacon*dt
+            ks = np.rint(shifts*f_beacon*dt)
 
             ks_per_loc[i] = ks
             maxima_per_loc[i] = maximum
 
         xx = np.array(xx)
         yy = np.array(yy)
-
-        best_idx = np.argmax(maxima_per_loc)
-        np.savetxt(__file__+f'.run{r}.txt')
-        print(ks_per_loc[best_idx])
-    
         if True: #plot maximum at test locations
             fig, axs = plt.subplots()
             axs.set_title(f"Grid Run {r}")
@@ -181,4 +183,31 @@ if __name__ == "__main__":
             fig.colorbar(sc, ax=axs)
 
             fig.savefig(__file__+f'.run{r}.pdf')
+
+        ## Save ks to file
+        best_idx = np.argmax(maxima_per_loc)
+        np.savetxt(__file__+f'.bestk.run{r}.txt', ks_per_loc[best_idx] )
+        print(ks_per_loc[best_idx])
+
+        ## Save maxima to file
+        np.savetxt(__file__+f'.maxima.run{r}.txt', maxima_per_loc)
+
+        # Abort if no improvement
+        if ( r!= 0 and (old_ks_per_loc == ks_per_loc[best_idx]).all() ):
+            print("No improvement, breaking")
+            break
+
+        # Prepare variables for next loop
+        old_ks_per_loc = ks_per_loc[best_idx]
+        xoff = xx[best_idx]
+        yoff = yy[best_idx]
+
+    # Save best ks to hdf5 antenna file
+    with h5py.File(antennas_fname, 'a') as fp:
+        group = fp.require_group('antennas')
+
+        for i, ant in enumerate(antennas):
+            h5ant = group[ant.name]
+            h5ant.attrs['best_k'] = old_ks_per_loc[i]
+
     plt.show()