Skip to content

Commit

Permalink
set efield position, fix efield reconstructor noise window definition
Browse files Browse the repository at this point in the history
  • Loading branch information
cg-laser committed Jun 5, 2024
1 parent cf1e900 commit 9a40073
Show file tree
Hide file tree
Showing 2 changed files with 19 additions and 12 deletions.
25 changes: 13 additions & 12 deletions NuRadioReco/modules/electricFieldSignalReconstructor.py
Original file line number Diff line number Diff line change
Expand Up @@ -61,16 +61,16 @@ def run(self, evt, station, det, debug=False):
electric_field[efp.signal_time] = signal_time

#
low_pos = int(130 * units.ns * electric_field.get_sampling_rate())
up_pos = int(210 * units.ns * electric_field.get_sampling_rate())
if(debug):
fig, ax = plt.subplots(1, 1)
sc = ax.scatter(trace_copy[1, low_pos:up_pos], trace_copy[2, low_pos:up_pos], c=electric_field.get_times()[low_pos:up_pos], s=5)
fig.colorbar(sc, ax=ax)
ax.set_aspect('equal')
ax.set_xlabel("eTheta")
ax.set_ylabel("ePhi")
fig.tight_layout()
# low_pos = int(130 * units.ns * electric_field.get_sampling_rate())
# up_pos = int(210 * units.ns * electric_field.get_sampling_rate())
# if(debug):
# fig, ax = plt.subplots(1, 1)
# sc = ax.scatter(trace_copy[1, low_pos:up_pos], trace_copy[2, low_pos:up_pos], c=electric_field.get_times()[low_pos:up_pos], s=5)
# fig.colorbar(sc, ax=ax)
# ax.set_aspect('equal')
# ax.set_xlabel("eTheta")
# ax.set_ylabel("ePhi")
# fig.tight_layout()

low_pos, up_pos = hp.get_interval(envelope_mag, scale=0.5)
v_start = trace_copy[:, signal_time_bin]
Expand Down Expand Up @@ -100,7 +100,8 @@ def run(self, evt, station, det, debug=False):
mask_signal_window = (times > (signal_time - self.__signal_window_pre)) & (times < (signal_time + self.__signal_window_post))
mask_noise_window = np.zeros_like(mask_signal_window, dtype=bool)
if(self.__noise_window > 0):
mask_noise_window[int(np.round((-self.__noise_window - 141.) * electric_field.get_sampling_rate())):int(np.round(-141. * electric_field.get_sampling_rate()))] = np.ones(int(np.round(self.__noise_window * electric_field.get_sampling_rate())), dtype=bool) # the last n bins
# set the noise window to the first "self.__noise_window" ns of the trace. If this cuts into the signal window, the noise window is reduced to not overlap with the signal window
mask_noise_window = times < min(times[0] + self.__noise_window, signal_time - self.__signal_window_pre)

signal_energy_fluence = trace_utilities.get_electric_field_energy_fluence(trace, times, mask_signal_window, mask_noise_window)
dt = times[1] - times[0]
Expand All @@ -127,7 +128,7 @@ def run(self, evt, station, det, debug=False):

# compute expeted polarization
if det is not None:
site = det.get_site(station.get_id())
site = det.get_site(station.get_id()).lower()
exp_efield = hp.get_lorentzforce_vector(electric_field[efp.zenith], electric_field[efp.azimuth], hp.get_magnetic_field_vector(site))
cs = coordinatesystems.cstrafo(electric_field[efp.zenith], electric_field[efp.azimuth], site=site)
exp_efield_onsky = cs.transform_from_ground_to_onsky(exp_efield)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,11 @@ def run(self, evt, station, det):

group_ids = select_channels_per_station(det, station.get_id(), station.get_channel_ids())
for gid, use_channels in group_ids.items():
pos = []
for channel_id in use_channels:
pos.append(det.get_relative_position(station.get_id(), channel_id))
pos = np.array(pos)
pos = np.average(pos, axis=0)
efield_antenna_factor = trace_utilities.get_efield_antenna_factor(station, frequencies, use_channels, det,
zenith, azimuth, self.antenna_provider)
V = np.zeros((len(use_channels), len(frequencies)), dtype=complex)
Expand All @@ -90,6 +95,7 @@ def run(self, evt, station, det):
E2[mask] = (V[-1] - efield_antenna_factor[-1][0] * E1)[mask] / efield_antenna_factor[-1][1][mask]

efield = ef.ElectricField(use_channels)
efield.set_position(pos)
efield.set_frequency_spectrum(np.array([np.zeros_like(E1), E1, E2]), sampling_rate)

efield.set_trace_start_time(station.get_channel(use_channels[0]).get_trace_start_time())
Expand Down

0 comments on commit 9a40073

Please sign in to comment.