ZH: ba measure phase: mask strong peak for every trace

simulations/airshower_beacon_simulation/ba_measure_beacon_phase.py

@@ -23,21 +23,28 @@ if __name__ == "__main__":
 
     from scriptlib import MyArgumentParser
     parser = MyArgumentParser()
-    parser.add_argument('--no-use-AxB', dest='use_AxB_trace', action='store_false')
+
+    group1 = parser.add_mutually_exclusive_group()
+    group1.add_argument('--AxB', dest='use_AxB_trace', action='store_true', help='Only use AxB trace, if both AxB and beacon are not used, we use the antenna polarisations.')
+    group1.add_argument('--beacon', dest='use_beacon_trace', action='store_true', help='Only use the beacon trace')
+
+    parser.add_argument('--N-mask', type=float, default=500, help='Mask N_MASK samples around the absolute maximum of the trace. (Default: %(default)d)')
+
     args = parser.parse_args()
 
     f_beacon_band = (49e-3,55e-3) #GHz
     allow_frequency_fitting = False
     read_frequency_from_file = True
+    N_mask = int(args.N_mask)
 
-    use_AxB_trace = args.use_AxB_trace
-    use_beacon_trace = True # only applicable if AxB = False
+    use_only_AxB_trace = args.use_AxB_trace
+    use_only_beacon_trace = args.use_beacon_trace # only applicable if AxB = False
 
     show_plots = args.show_plots
 
     figsize = (12,8)
 
-    print("use_AxB_trace:", use_AxB_trace, "use_beacon_trace:",use_beacon_trace)
+    print("use_only_AxB_trace:", use_only_AxB_trace, "use_only_beacon_trace:", use_only_beacon_trace)
 
     ####
     fname_dir = args.data_dir
@@ -83,59 +90,64 @@ if __name__ == "__main__":
             h5ant = group[name]
 
             # use E_AxB only instead of polarisations
-            if use_AxB_trace:
-                if 'E_AxB' not in h5ant.keys():
-                    print(f"Antenna does not have 'E_AxB' in {name}")
+            if use_only_AxB_trace:
+                traces_key = 'E_AxB'
+                if traces_key not in h5ant.keys():
+                    print(f"Antenna does not have '{traces_key}' in {name}")
                     sys.exit(1)
 
-                traces = h5ant['E_AxB']
-
+                traces = h5ant[traces_key]
                 t_trace = traces[0]
                 test_traces = [ traces[1] ]
                 orients = ['E_AxB']
 
-                # TODO: refine masking
-                # use beacon but remove where E_AxB-Beacon != 0
-                if True:
-                    N_pre, N_post = 250, 250 # TODO: make this configurable
+            # Only beacon
+            elif use_only_beacon_trace:
+                traces_key = 'filtered_traces'
+                if traces_key not in h5ant.keys():
+                    print(f"Antenna file corrupted? no '{traces_key}' in {name}")
+                    sys.exit(1)
 
-                    max_idx = np.argmax(test_traces[0])
-
-                    low_idx = max(0, max_idx-N_pre)
-                    high_idx = min(len(t_trace), max_idx+N_post)
-
-                    t_mask = np.ones(len(t_trace), dtype=bool)
-                    t_mask[low_idx:high_idx] = False
-
-                    t_trace = t_trace[t_mask]
-                    for j, t in enumerate(test_traces):
-                        test_traces[j] = t[t_mask]
-                        orients[j] = orients[j] + ' masked'
+                traces = h5ant[traces_key]
+                t_trace = traces[0]
+                test_traces = [traces[4]]
+                orients = ['B']
 
             # use separate polarisations
             else:
-                if 'traces' not in h5ant.keys():
-                    print(f"Antenna file corrupted? no 'traces' in {name}")
+                traces_key = 'filtered_traces'
+                if traces_key not in h5ant.keys():
+                    print(f"Antenna file corrupted? no '{traces_key}' in {name}")
                     sys.exit(1)
 
-                traces = h5ant['traces']
+                traces = h5ant[traces_key]
                 t_trace = traces[0]
+                test_traces = [traces[i] for i in range(1,4)]
+                orients = ['Ex', 'Ey', 'Ez']
 
-                if use_beacon_trace:
-                    # only take the Beacon trace
-                    test_traces = [traces[4]]
-                    orients = ['B']
-                else:
-                    test_traces = traces[1:]
-                    orients = ['Ex', 'Ey', 'Ez', 'B']
-
-            # modify the length of the traces
+            # Really only select the first component
             if False:
-                t_trace = t_trace[:len(t_trace)//2]
-                half_traces = []
-                for trace in test_traces:
-                    half_traces.append( trace[:len(trace)//2])
-                test_traces = half_traces
+                test_traces = [test_traces[0]]
+                orients = [orients[0]]
+
+            # TODO: refine masking
+            # use beacon but remove where E_AxB-Beacon != 0
+            # Uses the first traces as reference
+            if N_mask and orients[0] != 'B':
+                N_pre, N_post = N_mask//2, N_mask//2
+
+                max_idx = np.argmax(test_traces[0])
+
+                low_idx = max(0, max_idx-N_pre)
+                high_idx = min(len(t_trace), max_idx+N_post)
+
+                t_mask = np.ones(len(t_trace), dtype=bool)
+                t_mask[low_idx:high_idx] = False
+
+                t_trace = t_trace[t_mask]
+                for j, t in enumerate(test_traces):
+                    test_traces[j] = t[t_mask]
+                    orients[j] = orients[j] + ' masked'
 
             # Do Fourier Transforms
             # to find phases and amplitudes
@@ -155,7 +167,7 @@ if __name__ == "__main__":
                 amps = [ 3e-7 ]
 
             # choose highest amp
-            idx = 0
+            idx = np.argmax(amps)
             if False and len(beacon_phases) > 1:
                 #idx = np.argmax(amplitudes, axis=-1)
                 raise NotImplementedError