Skip to content

Commit

Permalink
Merge pull request #590 from nu-radio/hotfix/no-channel-0
Browse files Browse the repository at this point in the history 
fix some modules assuming the first channel has id 0
  • Loading branch information
@sjoerd-bouma
sjoerd-bouma authored Nov 16, 2023
2 parents e642c8b + a6bf9f2 commit 1504400
Show file tree
Hide file tree
Showing 6 changed files with 35 additions and 34 deletions.
4 changes: 2 additions & 2 deletions NuRadioReco/modules/beamFormingDirectionFitter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ def get_array_of_channels(station, det, zenith, azimuth, polarization):
time_shifts = np.zeros(8)
t_geos = np.zeros(8)

sampling_rate = station.get_channel(0).get_sampling_rate()
sampling_rate = next(station.iter_channels()).get_sampling_rate()
station_id = station.get_id()
site = det.get_site(station_id)
for iCh, channel in enumerate(station.get_electric_fields()):
Expand Down Expand Up @@ -185,7 +185,7 @@ def ll_regular_station(angles, evt, station, det, polarization, sampling_rate, p
positions = []
for chan in channels:
positions.append(det.get_relative_position(station_id, chan))
sampling_rate = station.get_channel(0).get_sampling_rate()
sampling_rate = station.get_channel(channels[0]).get_sampling_rate()

ll = opt.brute(
ll_regular_station,
Expand Down
2 changes: 1 addition & 1 deletion NuRadioReco/modules/correlationDirectionFitter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ def ll_regular_station_fft(angles, corr_02_fft, corr_13_fft, sampling_rate, posi
station_id = station.get_id()
positions_pairs = [[det.get_relative_position(station_id, channel_pairs[0][0]), det.get_relative_position(station_id, channel_pairs[0][1])],
[det.get_relative_position(station_id, channel_pairs[1][0]), det.get_relative_position(station_id, channel_pairs[1][1])]]
sampling_rate = station.get_channel(0).get_sampling_rate() # assume that channels have the same sampling rate
sampling_rate = station.get_channel(channel_pairs[0][0]).get_sampling_rate() # assume that channels have the same sampling rate
trace_start_time_pairs = [[station.get_channel(channel_pairs[0][0]).get_trace_start_time(), station.get_channel(channel_pairs[0][1]).get_trace_start_time()],
[station.get_channel(channel_pairs[1][0]).get_trace_start_time(), station.get_channel(channel_pairs[1][1]).get_trace_start_time()]]
# determine automatically if one channel has an inverted waveform with respect to the other
Expand Down
40 changes: 20 additions & 20 deletions NuRadioReco/modules/sphericalWaveFitter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,25 +13,25 @@

class sphericalWaveFitter:
" Fits position x,y, z of a source using spherical fit to channels "

def __init__(self):
pass


def begin(self, channel_ids = [0, 3, 9, 10]):
self.__channel_ids = channel_ids
pass


def run(self, evt, station, det, start_pulser_position, n_index = None, debug = True):

print("channels used for this reconstruction:", self.__channel_ids)


def get_distance(x, y):
return np.sqrt((x[0] - y[0])**2+ (x[1] - y[1])**2 + (x[2] - y[2])**2)


def get_time_delay_spherical_wave(position, ch_pair, n=n_index):
T0 = get_distance(position, det.get_relative_position(station_id, ch_pair[0]))/(constants.c/n_index)*units.s
T1 = get_distance(position , det.get_relative_position(station_id, ch_pair[1]))/(constants.c/n_index)*units.s
Expand All @@ -45,11 +45,11 @@ def get_time_delay_spherical_wave(position, ch_pair, n=n_index):
for j in range(i + 1, len(self.__channel_ids)):
relative_positions = det.get_relative_position(station_id, self.__channel_ids[i]) - det.get_relative_position(station_id, self.__channel_ids[j])
self.__relative_positions.append(relative_positions)

self.__channel_pairs.append([self.__channel_ids[i], self.__channel_ids[j]])
self.__sampling_rate = station.get_channel(0).get_sampling_rate()


self.__sampling_rate = station.get_channel(self.__channel_ids[0]).get_sampling_rate()
if debug:
fig, ax = plt.subplots( len(self.__channel_pairs), 2)

Expand All @@ -73,24 +73,24 @@ def likelihood(pulser_position, x = None, y = None, debug_corr = False):
ax.set_ylim((0, max(self.__correlation[ich])))
ax.axvline(pos, label = 'reconstruction', lw = 1, color = 'orange')
ax.axvline(self._pos_starting[ich], label = 'starting pos', lw = 1, color = 'green')
ax.set_title("channel pair {}".format( ch_pair), fontsize = 5)
ax.set_title("channel pair {}".format( ch_pair), fontsize = 5)
ax.legend(fontsize = 5)

if debug_corr:
fig.tight_layout()
fig.savefig("debug.pdf")

return -1*corr



trace = np.copy(station.get_channel(self.__channel_pairs[0][0]).get_trace())
self.__correlation = np.zeros((len(self.__channel_pairs), len(np.abs(scipy.signal.correlate(trace, trace))) ))

for ich, ch_pair in enumerate(self.__channel_pairs):
trace1 = np.copy(station.get_channel(self.__channel_pairs[ich][0]).get_trace())
trace2 = np.copy(station.get_channel(self.__channel_pairs[ich][1]).get_trace())

t_max1 = station.get_channel(self.__channel_pairs[ich][0]).get_times()[np.argmax(np.abs(trace1))]
t_max2 = station.get_channel(self.__channel_pairs[ich][1]).get_times()[np.argmax(np.abs(trace2))]
corr_range = 50 * units.ns
Expand All @@ -101,7 +101,7 @@ def likelihood(pulser_position, x = None, y = None, debug_corr = False):
else:
trace2[np.abs(station.get_channel(self.__channel_pairs[ich][1]).get_times() - t_max2) > corr_range] = 0
self.__correlation[ich] = np.abs(scipy.signal.correlate(trace1, trace2))



#### set positions for starting position ####
Expand All @@ -124,7 +124,7 @@ def likelihood(pulser_position, x = None, y = None, debug_corr = False):
print("reconstructed position: {}".format([ll[0], ll[1], ll[2]]))

if debug:

method = 'Nelder-Mead'
x = np.arange(x_start -3, x_start +3, dx)
y = np.arange(y_start - 3, y_start +3, dy)
Expand All @@ -137,7 +137,7 @@ def likelihood(pulser_position, x = None, y = None, debug_corr = False):
c = opt.minimize(likelihood, x0 = (start_pulser_position[2]), args = (x_i, y_i, False), method = method)
zz[ix, iy] = c.fun
zz_values[ix, iy] = c.x[0]

fig = plt.figure(figsize = (10, 5))
ax1 = fig.add_subplot(121)
pax1 = ax1.pcolor(xx, yy, zz.T)
Expand Down Expand Up @@ -170,4 +170,4 @@ def likelihood(pulser_position, x = None, y = None, debug_corr = False):

def end(self):
pass

9 changes: 5 additions & 4 deletions NuRadioReco/modules/voltageToAnalyticEfieldConverter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -314,10 +314,10 @@ def run(self, evt, station, det, debug=False, debug_plotpath=None,
efield_antenna_factor, V, V_timedomain = get_array_of_channels(station, use_channels,
det, zenith, azimuth, self.antenna_provider,
time_domain=True)
sampling_rate = station.get_channel(0).get_sampling_rate()
sampling_rate = station.get_channel(use_channels[0]).get_sampling_rate()
n_samples_time = V_timedomain.shape[1]

noise_RMS = det.get_noise_RMS(station.get_id(), 0)
noise_RMS = det.get_noise_RMS(station.get_id(), use_channels[0])

def obj_xcorr(params):
if(len(params) == 3):
Expand Down Expand Up @@ -482,9 +482,10 @@ def obj_amplitude_second_order(params, slope, phase, pos, compare='hilbert', deb
ax[iCh][1].axvline(imin, linestyle='--', alpha=.8)
ax[iCh][1].axvline(imax, linestyle='--', alpha=.8)
if(debug_obj and self.i_slope_fit_iterations % 50 == 0):
sim_channel = station.get_sim_station().get_channel(0)[0]
sim_station = station.get_sim_station()
sim_channel = next(sim_station.iter_channels())
ax[4][0].plot(sim_channel.get_frequencies() / units.MHz, np.abs(pulse.get_analytic_pulse_freq(ampTheta, slope, phase, len(sim_channel.get_times()), sim_channel.get_sampling_rate(), bandpass=bandpass, quadratic_term=second_order)), '--', color='orange')
ax[4][0].plot(sim_channel.get_frequencies() / units.MHz, np.abs(station.get_sim_station().get_channel(0)[0].get_frequency_spectrum()[1]), color='blue')
ax[4][0].plot(sim_channel.get_frequencies() / units.MHz, np.abs(sim_channel.get_frequency_spectrum()[1]), color='blue')
ax[4][1].plot(sim_channel.get_frequencies() / units.MHz, np.abs(pulse.get_analytic_pulse_freq(ampPhi, slope, phase, len(sim_channel.get_times()), sim_channel.get_sampling_rate(), bandpass=bandpass, quadratic_term=second_order)), '--', color='orange')
ax[4][1].plot(sim_channel.get_frequencies() / units.MHz, np.abs(sim_channel.get_frequency_spectrum()[2]), color='blue')
ax[4][0].set_xlim([20, 500])
Expand Down
10 changes: 5 additions & 5 deletions NuRadioReco/modules/voltageToEfieldConverter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ def get_array_of_channels(station, use_channels, det, zenith, azimuth,
tmax = time_shifts.max()
logger.debug("adding relative station time = {:.0f}ns".format((t_cables.min() + t_geos.max()) / units.ns))
logger.debug("delta t is {:.2f}".format(delta_t / units.ns))
trace_length = station.get_channel(0).get_times()[-1] - station.get_channel(0).get_times()[0]
trace_length = station.get_channel(use_channels[0]).get_times()[-1] - station.get_channel(use_channels[0]).get_times()[0]
debug_cut = 0
if(debug_cut):
fig, ax = plt.subplots(len(use_channels), 1)
Expand Down Expand Up @@ -110,11 +110,11 @@ def stacked_lstsq(L, b, rcond=1e-10):

class voltageToEfieldConverter:
"""
This module reconstructs the electric field by solving the system of equation that related the incident electric field via the antenna response functions to the measured voltages
This module reconstructs the electric field by solving the system of equation that related the incident electric field via the antenna response functions to the measured voltages
(see Eq. 4 of the NuRadioReco paper https://link.springer.com/article/10.1140/epjc/s10052-019-6971-5).
The module assumed that the electric field is identical at the antennas/channels that are used for the reconstruction. Furthermore, at least two antennas with
orthogonal polarization response are needed to reconstruct the 3dim electric field.
Alternatively, the polarization of the resulting efield could be forced to a single polarization component. In that case, a single antenna is sufficient.
orthogonal polarization response are needed to reconstruct the 3dim electric field.
Alternatively, the polarization of the resulting efield could be forced to a single polarization component. In that case, a single antenna is sufficient.
"""

def __init__(self):
Expand Down Expand Up @@ -183,7 +183,7 @@ def run(self, evt, station, det, use_channels=None, use_MC_direction=False, forc
efield3_f[1]])

electric_field = NuRadioReco.framework.electric_field.ElectricField(use_channels, [0, 0, 0])
electric_field.set_frequency_spectrum(efield3_f, station.get_channel(0).get_sampling_rate())
electric_field.set_frequency_spectrum(efield3_f, station.get_channel(use_channels[0]).get_sampling_rate())
electric_field.set_parameter(efp.zenith, zenith)
electric_field.set_parameter(efp.azimuth, azimuth)
# figure out the timing of the E-field
Expand Down
4 changes: 2 additions & 2 deletions NuRadioReco/modules/voltageToEfieldConverterPerChannel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,12 +53,12 @@ def run(self, evt, station, det, pol=0):
zenith = station[stnp.zenith]
azimuth = station[stnp.azimuth]

frequencies = station.get_channel(0).get_frequencies() # assuming that all channels have the same sampling rate and length
use_channels = det.get_channel_ids(station.get_id())
frequencies = station.get_channel(use_channels[0]).get_frequencies() # assuming that all channels have the same sampling rate and length
efield_antenna_factor = trace_utilities.get_efield_antenna_factor(station, frequencies, use_channels, det,
zenith, azimuth, self.antenna_provider)

sampling_rate = station.get_channel(0).get_sampling_rate()
sampling_rate = station.get_channel(use_channels[0]).get_sampling_rate()

for iCh, channel in enumerate(station.iter_channels()):
efield = ef.ElectricField([iCh])
Expand Down

0 comments on commit 1504400

Please sign in to comment.