diff --git a/NuRadioReco/modules/beamFormingDirectionFitter.py b/NuRadioReco/modules/beamFormingDirectionFitter.py index 5cc431d0a..126596baa 100644 --- a/NuRadioReco/modules/beamFormingDirectionFitter.py +++ b/NuRadioReco/modules/beamFormingDirectionFitter.py @@ -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()): @@ -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, diff --git a/NuRadioReco/modules/correlationDirectionFitter.py b/NuRadioReco/modules/correlationDirectionFitter.py index 5c6ca4d34..9cb2d8134 100644 --- a/NuRadioReco/modules/correlationDirectionFitter.py +++ b/NuRadioReco/modules/correlationDirectionFitter.py @@ -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 diff --git a/NuRadioReco/modules/sphericalWaveFitter.py b/NuRadioReco/modules/sphericalWaveFitter.py index ce9b8b943..a08c2821d 100644 --- a/NuRadioReco/modules/sphericalWaveFitter.py +++ b/NuRadioReco/modules/sphericalWaveFitter.py @@ -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 @@ -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) @@ -73,16 +73,16 @@ 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))) )) @@ -90,7 +90,7 @@ def likelihood(pulser_position, x = None, y = None, debug_corr = False): 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 @@ -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 #### @@ -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) @@ -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) @@ -170,4 +170,4 @@ def likelihood(pulser_position, x = None, y = None, debug_corr = False): def end(self): pass - + diff --git a/NuRadioReco/modules/voltageToAnalyticEfieldConverter.py b/NuRadioReco/modules/voltageToAnalyticEfieldConverter.py index d517b8d22..33bd3cff5 100644 --- a/NuRadioReco/modules/voltageToAnalyticEfieldConverter.py +++ b/NuRadioReco/modules/voltageToAnalyticEfieldConverter.py @@ -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): @@ -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]) diff --git a/NuRadioReco/modules/voltageToEfieldConverter.py b/NuRadioReco/modules/voltageToEfieldConverter.py index df6144a48..5c2f5fb08 100644 --- a/NuRadioReco/modules/voltageToEfieldConverter.py +++ b/NuRadioReco/modules/voltageToEfieldConverter.py @@ -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) @@ -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): @@ -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 diff --git a/NuRadioReco/modules/voltageToEfieldConverterPerChannel.py b/NuRadioReco/modules/voltageToEfieldConverterPerChannel.py index bdc79caf8..f55f5f875 100644 --- a/NuRadioReco/modules/voltageToEfieldConverterPerChannel.py +++ b/NuRadioReco/modules/voltageToEfieldConverterPerChannel.py @@ -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])