Merge commit 'remove-beacon-before-E_AxB' into main

simulations/airshower_beacon_simulation/ca_period_from_shower.py

@@ -189,6 +189,7 @@ if __name__ == "__main__":
     fname_dir = args.data_dir
     antennas_fname = path.join(fname_dir, beacon.antennas_fname)
     time_diffs_fname = 'time_diffs.hdf5' if not True else antennas_fname
+    tx_fname = path.join(fname_dir, beacon.tx_fname)
 
     ## This is a file indicating whether the k-finding algorithm was
     ## stopped early. This happens when the ks do not change between
@@ -203,6 +204,7 @@ if __name__ == "__main__":
 
     # Read in antennas from file
     _, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
+    _, __, txdata = beacon.read_tx_file(tx_fname)
     # Read original REvent
     ev = REvent(args.input_fname)
     # .. patch in our antennas
@@ -216,53 +218,81 @@ if __name__ == "__main__":
     freq_name = next(iter(freq_names))
     f_beacon = ev.antennas[0].beacon_info[freq_name]['freq']
 
-    # Prepare polarisation and passbands
-    rit.set_pol_and_bp(ev, low=low_bp, high=high_bp)
-
     ##
     ## Manipulate time and traces of each antenna
     ##
     ### Remove time due to true phase
     ### and optionally remove the beacon
+
+    ### Note: there is no use in changing *_AxB variables here (except for plotting),
+    ### they're recomputed by the upcoming rit.set_pol_and_bp call.
     measured_repair_offsets = beacon.read_antenna_clock_repair_offsets(ev.antennas, mode='phases', freq_name=freq_name)
     for i, ant in enumerate(ev.antennas):
-        ev.antennas[i].orig_t = ev.antennas[i].t_AxB
+        ev.antennas[i].orig_t = ev.antennas[i].t
+        ev.antennas[i].t += measured_repair_offsets[i]
+        # t_AxB will be set by the rit.set_pol_and_bp function
         ev.antennas[i].t_AxB += measured_repair_offsets[i]
 
         if apply_signal_window_from_max:
             N_pre, N_post = 250, 250 # TODO: make this configurable
 
-            max_idx = np.argmax(ant.E_AxB)
+            # Get max idx from all the traces
+            # and select the strongest
+            max_idx = []
+            maxs = []
 
+            for trace in [ant.Ex, ant.Ey, ant.Ez]:
+                idx = np.argmax(np.abs(trace))
+                max_idx.append(idx)
+                maxs.append( np.abs(trace[idx]) )
+
+            idx = np.argmax(maxs)
+            max_idx = max_idx[idx]
+
+            # Create window around max_idx
             low_idx = max(0, max_idx-N_pre)
             high_idx = min(len(ant.t), max_idx+N_post)
 
             ev.antennas[i].orig_t = ant.orig_t[low_idx:high_idx]
             ev.antennas[i].t = ant.t[low_idx:high_idx]
-            ev.antennas[i].t_AxB = ant.t_AxB[low_idx:high_idx]
 
             ev.antennas[i].Ex = ant.Ex[low_idx:high_idx]
             ev.antennas[i].Ey = ant.Ey[low_idx:high_idx]
             ev.antennas[i].Ez = ant.Ez[low_idx:high_idx]
+
+            ev.antennas[i].t_AxB = ant.t_AxB[low_idx:high_idx]
             ev.antennas[i].E_AxB = ant.E_AxB[low_idx:high_idx]
 
-
+        # .. and remove the beacon from the traces
+        # Note: ant.E_AxB is recalculated by rit.set_pol_and_bp
         if remove_beacon_from_trace:
             clock_phase = measured_repair_offsets[i]*2*np.pi*f_beacon
 
             beacon_phase = ant.beacon_info[freq_name]['beacon_phase']
             f = ant.beacon_info[freq_name]['freq']
             ampl = ant.beacon_info[freq_name]['amplitude']
-            calc_beacon = lib.sine_beacon(f, ev.antennas[i].t_AxB, amplitude=ampl, phase=beacon_phase-clock_phase)
+            calc_beacon = lib.sine_beacon(f, ev.antennas[i].t, amplitude=ampl, phase=beacon_phase-clock_phase)
+            tx_amps = txdata['amplitudes']
+            tx_amps_sum = np.sum(tx_amps)
 
+            # Split up contribution to the various polarisations
+            for j, amp in enumerate(tx_amps):
+                if j == 0:
+                    ev.antennas[i].Ex -= amp*(1/tx_amps_sum)*calc_beacon
+                elif j == 1:
+                    ev.antennas[i].Ey -= amp*(1/tx_amps_sum)*calc_beacon
+                elif j == 2:
+                    ev.antennas[i].Ez -= amp*(1/tx_amps_sum)*calc_beacon
+
+            # Subtract the beacon from E_AxB
             ev.antennas[i].E_AxB -= calc_beacon
 
             # Make a figure of the manipulated traces
-            if i == 2:
+            if i == 72:
                 orig_beacon_amplifier = ampl/max(ant.beacon)
 
                 fig, ax = plt.subplots(figsize=figsize)
-                ax.set_title(f"Signal and Beacon traces Antenna {i}")
+                ax.set_title(f"Signal and Beacon traces Antenna {ant.name}")
                 ax.set_xlabel("Time [ns]")
                 ax.set_ylabel("Amplitude [$\\mu V/m$]")
 
@@ -280,25 +310,27 @@ if __name__ == "__main__":
                         old_xlim = ax.get_xlim()
 
                         if True: # zoomed on part without peak of this trace
-                            wx, x = 100, 0#ant.t_AxB[np.argmax(ant.E_AxB)]
+                            wx, x = 200, min(ant.t_AxB)
-                            ax.set_xlim(x-wx, x+wx)
+                            ax.set_xlim(x-5, x+wx)
 
-                            fig.savefig(path.join(fig_dir, __file__+f'.traces.A{i}.zoomed.beacon.pdf'))
+                            fig.savefig(path.join(fig_dir, __file__+f'.traces.A{ant.name}.zoomed.beacon.pdf'))
 
                         if True: # zoomed on peak of this trace
                             idx = np.argmax(ev.antennas[i].E_AxB)
                             x = ev.antennas[i].t_AxB[idx]
-                            wx = 100
+                            wx = 300
-                            ax.set_xlim(x-wx, x+wx)
+                            ax.set_xlim(x-wx//2, x+wx//2)
-                            fig.savefig(path.join(fig_dir, __file__+f".traces.A{i}.zoomed.peak.pdf"))
+                            fig.savefig(path.join(fig_dir, __file__+f".traces.A{ant.name}.zoomed.peak.pdf"))
 
                         ax.set_xlim(*old_xlim)
 
-                    fig.savefig(path.join(fig_dir, __file__+f'.traces.A{i}.pdf'))
+                    fig.savefig(path.join(fig_dir, __file__+f'.traces.A{ant.name}.pdf'))
 
     if show_plots:
         plt.show()
 
+    # Prepare polarisation and passbands
+    rit.set_pol_and_bp(ev, low=low_bp, high=high_bp)
 
     # determine allowable ks per location
     dt = ev.antennas[0].t_AxB[1] - ev.antennas[0].t_AxB[0]

simulations/airshower_beacon_simulation/da_reconstruction.py

@@ -83,6 +83,36 @@ if __name__ == "__main__":
         ev.antennas[i].t += measured_repair_offsets[i]
         ev.antennas[i].t_AxB += measured_repair_offsets[i]
 
+        if apply_signal_window_from_max:
+            N_pre, N_post = 250, 250 # TODO: make this configurable
+
+            # Get max idx from all the traces
+            # and select the strongest
+            max_idx = []
+            maxs = []
+
+            for trace in [ant.Ex, ant.Ey, ant.Ez]:
+                idx = np.argmax(np.abs(trace))
+                max_idx.append(idx)
+                maxs.append( np.abs(trace[idx]) )
+
+            idx = np.argmax(maxs)
+            max_idx = max_idx[idx]
+
+            # Create window around max_idx
+            low_idx = max(0, max_idx-N_pre)
+            high_idx = min(len(ant.t), max_idx+N_post)
+
+            ev.antennas[i].orig_t = ant.orig_t[low_idx:high_idx]
+            ev.antennas[i].t = ant.t[low_idx:high_idx]
+
+            ev.antennas[i].Ex = ant.Ex[low_idx:high_idx]
+            ev.antennas[i].Ey = ant.Ey[low_idx:high_idx]
+            ev.antennas[i].Ez = ant.Ez[low_idx:high_idx]
+
+            ev.antennas[i].t_AxB = ant.t_AxB[low_idx:high_idx]
+            ev.antennas[i].E_AxB = ant.E_AxB[low_idx:high_idx]
+
         # .. and remove the beacon from the traces
         # Note: ant.E_AxB is recalculated by rit.set_pol_and_bp
         if remove_beacon_from_traces:
@@ -111,7 +141,7 @@ if __name__ == "__main__":
             ##
             ## Make a figure of the manipulated traces
             ##
-            if i == 2:
+            if i == 72:
                 orig_beacon_amplifier = ampl_AxB/max(ant.beacon)
 
                 for k in range(2):
@@ -127,7 +157,7 @@ if __name__ == "__main__":
                         fname_extra = ".Ex"
 
                     fig, ax = plt.subplots(figsize=figsize)
-                    ax.set_title(f"Signal and Beacon traces Antenna {i}")
+                    ax.set_title(f"Signal and Beacon traces Antenna {ant.name}")
                     ax.set_xlabel("Time [ns]")
                     ax.set_ylabel("Amplitude [$\\mu V/m$]")
 
@@ -151,14 +181,14 @@ if __name__ == "__main__":
                                 wx, x = 100, 100
                                 ax.set_xlim(x-wx, x+wx)
 
-                                fig.savefig(path.join(fig_dir, __file__+f'.traces.A{i}.zoomed.beacon{fname_extra}.pdf'))
+                                fig.savefig(path.join(fig_dir, __file__+f'.traces.A{ant.name}.zoomed.beacon{fname_extra}.pdf'))
 
                             if True: # zoomed on peak of this trace
                                 idx = np.argmax(ev.antennas[i].E_AxB)
                                 x = ev.antennas[i].t_AxB[idx]
                                 wx = 100
                                 ax.set_xlim(x-wx, x+wx)
-                                fig.savefig(path.join(fig_dir, __file__+f".traces.A{i}.zoomed.peak{fname_extra}.pdf"))
+                                fig.savefig(path.join(fig_dir, __file__+f".traces.A{ant.name}.zoomed.peak{fname_extra}.pdf"))
 
                             ax.set_xlim(*old_xlim)
 

simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py

@@ -135,7 +135,7 @@ if __name__ == "__main__":
                 elif j == 2:
                     ev.antennas[i].Ez -= amp*(1/tx_amps_sum)*calc_beacon
 
-            # 
+            # Subtract the beacon from E_AxB
             ev.antennas[i].E_AxB -= calc_beacon
 
     # Slice the traces to a small part around the peak
@@ -143,7 +143,19 @@ if __name__ == "__main__":
         N_pre, N_post = 250, 250 # TODO: make this configurable
 
         for i, ant in enumerate(ev.antennas):
-            max_idx = np.argmax(ant.E_AxB)
+            # Get max idx from all the traces
+            # and select the strongest
+            max_idx = []
+            maxs = []
+
+            for trace in [ant.Ex, ant.Ey, ant.Ez]:
+                idx = np.argmax(np.abs(trace))
+                max_idx.append(idx)
+                maxs.append( np.abs(trace[idx]) )
+
+            idx = np.argmax(maxs)
+            max_idx = max_idx[idx]
+
 
             low_idx = max(0, max_idx-N_pre)
             high_idx = min(len(ant.t), max_idx+N_post)
@@ -156,13 +168,13 @@ if __name__ == "__main__":
             ev.antennas[i].Ez = ant.Ez[low_idx:high_idx]
             ev.antennas[i].E_AxB = ant.E_AxB[low_idx:high_idx]
 
+    ## Apply polarisation and bandpass filter
+    rit.set_pol_and_bp(ev)
+
     # backup antenna times
     backup_antenna_t = [ ant.t for ant in ev.antennas ]
     backup_antenna_t_AxB = [ ant.t_AxB for ant in ev.antennas ]
 
-    ## Apply polarisation and bandpass filter
-    rit.set_pol_and_bp(ev)
-
     with joblib.parallel_backend("loky"):
         for case in wanted_cases:
             print(f"Starting {case} figure")