diff --git a/NuRadioMC/EvtGen/generate_unforced.py b/NuRadioMC/EvtGen/generate_unforced.py index 2c228a68f..718f3fe0f 100644 --- a/NuRadioMC/EvtGen/generate_unforced.py +++ b/NuRadioMC/EvtGen/generate_unforced.py @@ -328,7 +328,7 @@ def points_in_cylinder(pt1, pt2, r, q): 'flavors': [], 'energies': []} # calculate rotation matrix to transform position on area to 3D - mask_int = np.zeros_like(mask, dtype=np.bool) + mask_int = np.zeros_like(mask, dtype=bool) t0 = time.perf_counter() n_cylinder = 0 for j, i in enumerate(np.arange(n_events, dtype=int)[mask]): @@ -586,7 +586,7 @@ def points_in_cylinder(pt1, pt2, r, q): data_sets['flavors'] = np.ones(np.sum(mask_int)) data_sets["event_ids"] = np.arange(np.sum(mask_int)) + start_event_id data_sets["n_interaction"] = np.ones(np.sum(mask_int), dtype=int) - data_sets["vertex_times"] = np.zeros(np.sum(mask_int), dtype=np.float) + data_sets["vertex_times"] = np.zeros(np.sum(mask_int), dtype=float) data_sets["interaction_type"] = inelasticities.get_ccnc(np.sum(mask_int)) data_sets["inelasticity"] = inelasticities.get_neutrino_inelasticity(np.sum(mask_int)) diff --git a/NuRadioMC/EvtGen/generator.py b/NuRadioMC/EvtGen/generator.py index d3a73ce65..ae68676dc 100644 --- a/NuRadioMC/EvtGen/generator.py +++ b/NuRadioMC/EvtGen/generator.py @@ -891,7 +891,7 @@ def generate_surface_muons(filename, n_events, Emin, Emax, data_sets["event_group_ids"] = np.arange(i_batch * max_n_events_batch, i_batch * max_n_events_batch + n_events_batch, dtype=int) + start_event_id data_sets["n_interaction"] = np.ones(n_events_batch, dtype=int) - data_sets["vertex_times"] = np.zeros(n_events_batch, dtype=np.float) + data_sets["vertex_times"] = np.zeros(n_events_batch, dtype=float) # generate neutrino flavors randomly @@ -928,7 +928,7 @@ def generate_surface_muons(filename, n_events, Emin, Emax, if('fiducial_rmax' in attributes): mask_phi = mask_arrival_azimuth(data_sets, attributes['fiducial_rmax']) # this currently only works for cylindrical volumes else: - mask_phi = np.ones(len(data_sets["event_group_ids"]), dtype=np.bool) + mask_phi = np.ones(len(data_sets["event_group_ids"]), dtype=bool) # TODO: combine with `get_intersection_volume_neutrino` function for iE, event_id in enumerate(data_sets["event_group_ids"]): if not mask_phi[iE]: @@ -1230,7 +1230,7 @@ def generate_eventlist_cylinder(filename, n_events, Emin, Emax, data_sets["event_group_ids"] = np.arange(i_batch * max_n_events_batch, i_batch * max_n_events_batch + n_events_batch) + start_event_id logger.debug("generating number of interactions") data_sets["n_interaction"] = np.ones(n_events_batch, dtype=int) - data_sets["vertex_times"] = np.zeros(n_events_batch, dtype=np.float) + data_sets["vertex_times"] = np.zeros(n_events_batch, dtype=float) # generate neutrino flavors randomly logger.debug("generating flavors") diff --git a/NuRadioMC/SignalGen/ARZ/ARZ.py b/NuRadioMC/SignalGen/ARZ/ARZ.py index 4b8523200..92827124c 100644 --- a/NuRadioMC/SignalGen/ARZ/ARZ.py +++ b/NuRadioMC/SignalGen/ARZ/ARZ.py @@ -63,7 +63,7 @@ def get_vector_potential( profile_ce: array of floats charge-excess values of the charge excess profile shower_type: string (default "HAD") - type of shower, either "HAD" (hadronic), "EM" (electromagnetic) or "TAU" (tau lepton induced) + type of shower, either "HAD" (hadronic) or "EM" (electromagnetic) n_index: float (default 1.78) index of refraction where the shower development takes place distance: float (default 1km) @@ -239,7 +239,7 @@ def get_dist_shower(X, z): v[0] = vperp_x v[1] = vperp_y v[2] = vperp_z -# v = np.array([vperp_x, vperp_y, vperp_z], dtype=np.float64) +# v = np.array([vperp_x, vperp_y, vperp_z], dtype=float) """ Function F_p Eq.(15) PRD paper. """ @@ -468,6 +468,37 @@ def set_interpolation_factor2(self, interp_factor): """ self._interp_factor2 = interp_factor + def get_shower_profile(self, shower_energy, shower_type, iN): + """ + Gets a charge-excess profile + + Parameters + ---------- + shower_energy: float + the energy of the shower + shower_type: string (default "HAD") + type of shower, either "HAD" (hadronic) or "EM" (electromagnetic) + iN: int + specify shower number + + Returns + ------- + depth, excess: two arrays of floats + slant depths and charge profile amplitudes + """ + + energies = np.array([*self._library[shower_type]]) + iE = np.argmin(np.abs(energies - shower_energy)) + + rescaling_factor = shower_energy / energies[iE] + + profiles = self._library[shower_type][energies[iE]] + profile_depth = profiles['depth'] + profile_ce = profiles['charge_excess'][iN] * rescaling_factor + + return profile_depth, profile_ce + + def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, shift_for_xmax=False, same_shower=False, iN=None, output_mode='trace', maximum_angle=20 * units.deg): """ @@ -483,12 +514,8 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, s number of samples in the time domain dt: float size of one time bin in units of time - profile_depth: array of floats - shower depth values of the charge excess profile - profile_ce: array of floats - charge-excess values of the charge excess profile shower_type: string (default "HAD") - type of shower, either "HAD" (hadronic), "EM" (electromagnetic) or "TAU" (tau lepton induced) + type of shower, either "HAD" (hadronic) or "EM" (electromagnetic) n_index: float (default 1.78) index of refraction where the shower development takes place R: float (default 1km) @@ -537,6 +564,7 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, s rescaling_factor = shower_energy / energies[iE] logger.info("shower energy of {:.3g}eV requested, closest available energy is {:.3g}eV. The amplitude of the charge-excess profile will be rescaled accordingly by a factor of {:.2f}".format(shower_energy / units.eV, energies[iE] / units.eV, rescaling_factor)) profiles = self._library[shower_type][energies[iE]] + N_profiles = len(profiles['charge_excess']) if(iN is None or np.isnan(iN)): @@ -597,7 +625,7 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, s logger.error("Tau showers are not yet implemented") raise NotImplementedError("Tau showers are not yet implemented") else: - msg = "showers of type {} are not implemented. Use 'HAD', 'EM' or 'TAU'".format(shower_type) + msg = "showers of type {} are not implemented. Use 'HAD', 'EM'".format(shower_type) logger.error(msg) raise NotImplementedError(msg) if self._use_numba: @@ -669,7 +697,7 @@ def get_vector_potential_fast(self, shower_energy, theta, N, dt, profile_depth, profile_ce: array of floats charge-excess values of the charge excess profile shower_type: string (default "HAD") - type of shower, either "HAD" (hadronic), "EM" (electromagnetic) or "TAU" (tau lepton induced) + type of shower, either "HAD" (hadronic) or "EM" (electromagnetic) n_index: float (default 1.78) index of refraction where the shower development takes place distance: float (default 1km) @@ -1136,7 +1164,7 @@ def get_time_trace(self, shower_energy, theta, N, dt, shower_type, n_index, R, dt: float size of one time bin in units of time shower_type: string (default "HAD") - type of shower, either "HAD" (hadronic), "EM" (electromagnetic) or "TAU" (tau lepton induced) + type of shower, either "HAD" (hadronic) or "EM" (electromagnetic) n_index: float (default 1.78) index of refraction where the shower development takes place R: float (default 1km) diff --git a/NuRadioMC/SignalProp/CPPAnalyticRayTracing/setup.py b/NuRadioMC/SignalProp/CPPAnalyticRayTracing/setup.py index 932437653..adc09a364 100644 --- a/NuRadioMC/SignalProp/CPPAnalyticRayTracing/setup.py +++ b/NuRadioMC/SignalProp/CPPAnalyticRayTracing/setup.py @@ -16,7 +16,7 @@ Extension('wrapper', ['wrapper.pyx'], include_dirs=[numpy.get_include(), '../../utilities/', str(os.environ['GSLDIR']) + '/include/'], library_dirs=[str(os.environ['GSLDIR']) + '/lib/'], - extra_compile_args=['-O3',"-mfpmath=sse"], + extra_compile_args=['-O3'], libraries=['gsl', 'gslcblas'], language='c++' ), diff --git a/NuRadioMC/SignalProp/analyticraytracing.py b/NuRadioMC/SignalProp/analyticraytracing.py index d9e3c670a..bac2c061b 100644 --- a/NuRadioMC/SignalProp/analyticraytracing.py +++ b/NuRadioMC/SignalProp/analyticraytracing.py @@ -639,7 +639,7 @@ def dt(t, C_0, frequency): mask = frequency > 0 freqs = self.__get_frequencies_for_attenuation(frequency, max_detector_freq) gamma_turn, z_turn = self.get_turning_point(self.medium.n_ice ** 2 - C_0 ** -2) - print("_use_optimized_calculation", self._use_optimized_calculation) + self.__logger.info("_use_optimized_calculation", self._use_optimized_calculation) if self._use_optimized_calculation: # The integration of the attenuation factor along the path with scipy.quad is inefficient. The @@ -1751,8 +1751,8 @@ def set_start_and_end_point(self, x1, x2): if(self._X2[2] < self._X1[2]): self._swap = True self.__logger.debug('swap = True') - self._X2 = np.array(x1, dtype =np.float) - self._X1 = np.array(x2, dtype =np.float) + self._X2 = np.array(x1, dtype =float) + self._X1 = np.array(x2, dtype =float) dX = self._X2 - self._X1 self._dPhi = -np.arctan2(dX[1], dX[0]) diff --git a/NuRadioMC/SignalProp/propagation_base_class.py b/NuRadioMC/SignalProp/propagation_base_class.py index 376f5d19d..0ba908394 100644 --- a/NuRadioMC/SignalProp/propagation_base_class.py +++ b/NuRadioMC/SignalProp/propagation_base_class.py @@ -92,8 +92,8 @@ def set_start_and_end_point(self, x1, x2): stop point of the ray """ self.reset_solutions() - self._X1 = np.array(x1, dtype =np.float) - self._X2 = np.array(x2, dtype = np.float) + self._X1 = np.array(x1, dtype =float) + self._X2 = np.array(x2, dtype = float) if (self._n_reflections): if (self._X1[2] < self._medium.reflection or self._X2[2] < self._medium.reflection): self.__logger.error("start or stop point is below the reflective bottom layer at {:.1f}m".format( diff --git a/NuRadioMC/examples/03_station_coincidences/A05visualization.py b/NuRadioMC/examples/03_station_coincidences/A05visualization.py index 99095c5c0..e69c6f2d5 100644 --- a/NuRadioMC/examples/03_station_coincidences/A05visualization.py +++ b/NuRadioMC/examples/03_station_coincidences/A05visualization.py @@ -79,7 +79,7 @@ x2 = det.get_relative_position(101, iC) if(plot): ax.plot([x2[0]], [x2[1]], [x2[2]], 'ko') - if(j != 0 and (~(np.array(triggered_deep[iE], dtype=np.bool) & mask))[j]): + if(j != 0 and (~(np.array(triggered_deep[iE], dtype=bool) & mask))[j]): continue vertex = np.array([fin['xx'][iE], fin['yy'][iE], fin['zz'][iE]]) # print(fin.keys()) diff --git a/NuRadioMC/simulation/scripts/T05visualize_sim_output.py b/NuRadioMC/simulation/scripts/T05visualize_sim_output.py index a1c046cbf..9585974db 100644 --- a/NuRadioMC/simulation/scripts/T05visualize_sim_output.py +++ b/NuRadioMC/simulation/scripts/T05visualize_sim_output.py @@ -42,14 +42,14 @@ plot_folder = os.path.join(dirname, 'plots', filename, args.trigger_name[0]) if(not os.path.exists(plot_folder)): os.makedirs(plot_folder) - triggered = np.zeros(len(fin['multiple_triggers'][:, 0]), dtype=np.bool) + triggered = np.zeros(len(fin['multiple_triggers'][:, 0]), dtype=bool) for trigger in args.trigger_name[1:]: iTrigger = np.squeeze(np.argwhere(fin.attrs['trigger_names'] == trigger)) - triggered = triggered | np.array(fin['multiple_triggers'][:, iTrigger], dtype=np.bool) + triggered = triggered | np.array(fin['multiple_triggers'][:, iTrigger], dtype=bool) else: trigger_name = args.trigger_name[0] iTrigger = np.argwhere(fin.attrs['trigger_names'] == trigger_name) - triggered = np.array(fin['multiple_triggers'][:, iTrigger], dtype=np.bool) + triggered = np.array(fin['multiple_triggers'][:, iTrigger], dtype=bool) print("\tyou selected '{}'".format(trigger_name)) plot_folder = os.path.join(dirname, 'plots', filename, trigger_name) if(not os.path.exists(plot_folder)): @@ -147,14 +147,14 @@ def a(theta): if(len(args.trigger_name) > 1): print("trigger {} selected which is a combination of {}".format(args.trigger_name[0], args.trigger_name[1:])) trigger_name = args.trigger_name[0] - triggered = np.zeros(len(station['multiple_triggers'][:, 0]), dtype=np.bool) + triggered = np.zeros(len(station['multiple_triggers'][:, 0]), dtype=bool) for trigger in args.trigger_name[1:]: iTrigger = np.squeeze(np.argwhere(fin.attrs['trigger_names'] == trigger)) - triggered = triggered | np.array(station['multiple_triggers'][:, iTrigger], dtype=np.bool) + triggered = triggered | np.array(station['multiple_triggers'][:, iTrigger], dtype=bool) else: trigger_name = args.trigger_name[0] iTrigger = np.argwhere(fin.attrs['trigger_names'] == trigger_name) - triggered = np.array(station['multiple_triggers'][:, iTrigger], dtype=np.bool) + triggered = np.array(station['multiple_triggers'][:, iTrigger], dtype=bool) print("\tyou selected '{}'".format(trigger_name)) ########################### diff --git a/NuRadioMC/simulation/simulation.py b/NuRadioMC/simulation/simulation.py index 35b6be0ac..42a1d19fa 100644 --- a/NuRadioMC/simulation/simulation.py +++ b/NuRadioMC/simulation/simulation.py @@ -306,7 +306,7 @@ def __init__(self, inputfilename, self._amplification_per_channel[self._station_id] = {} for channel_id in range(self._det.get_number_of_channels(self._station_id)): ff = np.linspace(0, 0.5 / self._dt, 10000) - filt = np.ones_like(ff, dtype=np.complex) + filt = np.ones_like(ff, dtype=complex) for i, (name, instance, kwargs) in enumerate(self._evt.iter_modules(self._station_id)): if hasattr(instance, "get_filter"): filt *= instance.get_filter(ff, self._station_id, channel_id, self._det, **kwargs) @@ -1150,12 +1150,12 @@ def _calculate_signal_properties(self): def _create_empty_multiple_triggers(self): if 'trigger_names' not in self._mout_attrs: self._mout_attrs['trigger_names'] = np.array([]) - self._mout['multiple_triggers'] = np.zeros((self._n_showers, 1), dtype=np.bool) + self._mout['multiple_triggers'] = np.zeros((self._n_showers, 1), dtype=bool) for station_id in self._station_ids: sg = self._mout_groups[station_id] n_showers = sg['launch_vectors'].shape[0] - sg['multiple_triggers'] = np.zeros((n_showers, 1), dtype=np.bool) - sg['triggered'] = np.zeros(n_showers, dtype=np.bool) + sg['multiple_triggers'] = np.zeros((n_showers, 1), dtype=bool) + sg['triggered'] = np.zeros(n_showers, dtype=bool) def _create_trigger_structures(self): @@ -1170,13 +1170,13 @@ def _create_trigger_structures(self): # simulated triggers is unknown at the beginning. So we check if the key already exists and if not, # we first create this data structure if 'multiple_triggers' not in self._mout: - self._mout['multiple_triggers'] = np.zeros((self._n_showers, len(self._mout_attrs['trigger_names'])), dtype=np.bool) + self._mout['multiple_triggers'] = np.zeros((self._n_showers, len(self._mout_attrs['trigger_names'])), dtype=bool) self._mout['trigger_times'] = np.nan * np.zeros_like(self._mout['multiple_triggers'], dtype=float) # for station_id in self._station_ids: # sg = self._mout_groups[station_id] -# sg['multiple_triggers'] = np.zeros((self._n_showers, len(self._mout_attrs['trigger_names'])), dtype=np.bool) +# sg['multiple_triggers'] = np.zeros((self._n_showers, len(self._mout_attrs['trigger_names'])), dtype=bool) elif extend_array: - tmp = np.zeros((self._n_showers, len(self._mout_attrs['trigger_names'])), dtype=np.bool) + tmp = np.zeros((self._n_showers, len(self._mout_attrs['trigger_names'])), dtype=bool) nx, ny = self._mout['multiple_triggers'].shape tmp[:, 0:ny] = self._mout['multiple_triggers'] self._mout['multiple_triggers'] = tmp @@ -1186,7 +1186,7 @@ def _create_trigger_structures(self): self._mout['trigger_times'] = tmp_t # for station_id in self._station_ids: # sg = self._mout_groups[station_id] -# tmp = np.zeros((self._n_showers, len(self._mout_attrs['trigger_names'])), dtype=np.bool) +# tmp = np.zeros((self._n_showers, len(self._mout_attrs['trigger_names'])), dtype=bool) # nx, ny = sg['multiple_triggers'].shape # tmp[:, 0:ny] = sg['multiple_triggers'] # sg['multiple_triggers'] = tmp @@ -1199,10 +1199,10 @@ def _save_triggers_to_hdf5(self, sg, local_shower_index, global_shower_index): # the information fo the current station and event group n_showers = sg['launch_vectors'].shape[0] if 'multiple_triggers' not in sg: - sg['multiple_triggers'] = np.zeros((n_showers, len(self._mout_attrs['trigger_names'])), dtype=np.bool) + sg['multiple_triggers'] = np.zeros((n_showers, len(self._mout_attrs['trigger_names'])), dtype=bool) sg['trigger_times'] = np.nan * np.zeros_like(sg['multiple_triggers'], dtype=float) elif extend_array: - tmp = np.zeros((n_showers, len(self._mout_attrs['trigger_names'])), dtype=np.bool) + tmp = np.zeros((n_showers, len(self._mout_attrs['trigger_names'])), dtype=bool) nx, ny = sg['multiple_triggers'].shape tmp[:, 0:ny] = sg['multiple_triggers'] sg['multiple_triggers'] = tmp @@ -1213,7 +1213,7 @@ def _save_triggers_to_hdf5(self, sg, local_shower_index, global_shower_index): self._output_event_group_ids[self._station_id].append(self._evt.get_run_number()) self._output_sub_event_ids[self._station_id].append(self._evt.get_id()) - multiple_triggers = np.zeros(len(self._mout_attrs['trigger_names']), dtype=np.bool) + multiple_triggers = np.zeros(len(self._mout_attrs['trigger_names']), dtype=bool) trigger_times = np.nan*np.zeros_like(multiple_triggers) for iT, trigger_name in enumerate(self._mout_attrs['trigger_names']): if self._station.has_trigger(trigger_name): @@ -1262,8 +1262,8 @@ def _create_meta_output_datastructures(self): self._mout = {} self._mout_attributes = {} self._mout['weights'] = np.zeros(self._n_showers) - self._mout['triggered'] = np.zeros(self._n_showers, dtype=np.bool) -# self._mout['multiple_triggers'] = np.zeros((self._n_showers, self._number_of_triggers), dtype=np.bool) + self._mout['triggered'] = np.zeros(self._n_showers, dtype=bool) +# self._mout['multiple_triggers'] = np.zeros((self._n_showers, self._number_of_triggers), dtype=bool) self._mout_attributes['trigger_names'] = None self._amplitudes = {} self._amplitudes_envelope = {} @@ -1289,7 +1289,7 @@ def _create_meta_output_datastructures(self): def _create_station_output_structure(self, n_showers, n_antennas): nS = self._raytracer.get_number_of_raytracing_solutions() # number of possible ray-tracing solutions sg = {} - sg['triggered'] = np.zeros(n_showers, dtype=np.bool) + sg['triggered'] = np.zeros(n_showers, dtype=bool) # we need the reference to the shower id to be able to find the correct shower in the upper level hdf5 file sg['shower_id'] = np.zeros(n_showers, dtype=int) * -1 sg['event_id_per_shower'] = np.zeros(n_showers, dtype=int) * -1 @@ -1424,7 +1424,7 @@ def _write_output_file(self, empty=False): # the multiple triggeres 2d array might have different number of entries per event # because the number of different triggers can increase dynamically # therefore we first create an array with the right size and then fill it - tmp = np.zeros((n_events_for_station, n_triggers), dtype=np.bool) + tmp = np.zeros((n_events_for_station, n_triggers), dtype=bool) for iE, values in enumerate(self._output_multiple_triggers_station[station_id]): tmp[iE] = values sg['multiple_triggers_per_event'] = tmp @@ -1452,10 +1452,10 @@ def _write_output_file(self, empty=False): fout["station_{:d}".format(station_id)].attrs['Vrms'] = list(self._Vrms_per_channel[station_id].values()) fout["station_{:d}".format(station_id)].attrs['bandwidth'] = list(self._bandwidth_per_channel[station_id].values()) - fout.attrs.create("Tnoise", self._noise_temp, dtype=np.float) - fout.attrs.create("Vrms", self._Vrms, dtype=np.float) - fout.attrs.create("dt", self._dt, dtype=np.float) - fout.attrs.create("bandwidth", self._bandwidth, dtype=np.float) + fout.attrs.create("Tnoise", self._noise_temp, dtype=float) + fout.attrs.create("Vrms", self._Vrms, dtype=float) + fout.attrs.create("dt", self._dt, dtype=float) + fout.attrs.create("bandwidth", self._bandwidth, dtype=float) fout.attrs['n_samples'] = self._n_samples fout.attrs['config'] = yaml.dump(self._cfg) diff --git a/NuRadioMC/test/SingleEvents/1e18_output_noise_reference.hdf5 b/NuRadioMC/test/SingleEvents/1e18_output_noise_reference.hdf5 index 093452baa..83cbb5536 100644 Binary files a/NuRadioMC/test/SingleEvents/1e18_output_noise_reference.hdf5 and b/NuRadioMC/test/SingleEvents/1e18_output_noise_reference.hdf5 differ diff --git a/NuRadioMC/test/SingleEvents/validate.sh b/NuRadioMC/test/SingleEvents/validate.sh index 5aa12b647..233784c04 100755 --- a/NuRadioMC/test/SingleEvents/validate.sh +++ b/NuRadioMC/test/SingleEvents/validate.sh @@ -1,9 +1,8 @@ #!/bin/bash -python T02RunSimulation.py 1e18_output_reference.hdf5 surface_station_1GHz.json config.yaml 1e18_output.hdf5 1e18_output.nur - -python T03validate.py 1e18_output.hdf5 1e18_output_reference.hdf5 +python3 T02RunSimulation.py 1e18_output_reference.hdf5 surface_station_1GHz.json config.yaml 1e18_output.hdf5 1e18_output.nur +python3 T03validate.py 1e18_output.hdf5 1e18_output_reference.hdf5 # cleanup -rm -v NuRadioMC/test/SingleEvents/1e18_output.hdf5 \ No newline at end of file +rm -v NuRadioMC/test/SingleEvents/1e18_output.hdf5 diff --git a/NuRadioMC/test/Veff/1e18eV/T03check_output_noise.py b/NuRadioMC/test/Veff/1e18eV/T03check_output_noise.py index a78fe4bfc..f1ebfc01f 100755 --- a/NuRadioMC/test/Veff/1e18eV/T03check_output_noise.py +++ b/NuRadioMC/test/Veff/1e18eV/T03check_output_noise.py @@ -11,7 +11,7 @@ # the event generation has a fixed seed and I switched to Alvarez2000 (also no randomness) # thus, the Veff has no statistical scatter -Veff_mean = 8.17491 +Veff_mean = 7.86364 Veff_sigma = 0.0001 path = os.path.dirname(os.path.abspath(__file__)) diff --git a/NuRadioMC/utilities/Veff.py b/NuRadioMC/utilities/Veff.py index 9ee5eaa14..00aecb933 100644 --- a/NuRadioMC/utilities/Veff.py +++ b/NuRadioMC/utilities/Veff.py @@ -335,35 +335,35 @@ def get_Veff_Aeff_single( return out for iT, trigger_name in enumerate(trigger_names): - triggered = np.array(fin['multiple_triggers'][:, iT], dtype=np.bool) + triggered = np.array(fin['multiple_triggers'][:, iT], dtype=bool) triggered = remove_duplicate_triggers(triggered, fin['event_group_ids']) out[veff_aeff][trigger_name] = get_veff_output(volume_proj_area, np.sum(weights[triggered]), n_events) for trigger_name, values in iteritems(trigger_combinations): indiv_triggers = values['triggers'] - triggered = np.zeros_like(fin['multiple_triggers'][:, 0], dtype=np.bool) + triggered = np.zeros_like(fin['multiple_triggers'][:, 0], dtype=bool) if isinstance(indiv_triggers, str): triggered = triggered | \ - np.array(fin['multiple_triggers'][:, trigger_names_dict[indiv_triggers]], dtype=np.bool) + np.array(fin['multiple_triggers'][:, trigger_names_dict[indiv_triggers]], dtype=bool) else: for indiv_trigger in indiv_triggers: triggered = triggered | \ - np.array(fin['multiple_triggers'][:, trigger_names_dict[indiv_trigger]], dtype=np.bool) + np.array(fin['multiple_triggers'][:, trigger_names_dict[indiv_trigger]], dtype=bool) if 'triggerAND' in values: triggered = triggered & \ - np.array(fin['multiple_triggers'][:, trigger_names_dict[values['triggerAND']]], dtype=np.bool) + np.array(fin['multiple_triggers'][:, trigger_names_dict[values['triggerAND']]], dtype=bool) if 'notriggers' in values: indiv_triggers = values['notriggers'] if(isinstance(indiv_triggers, str)): triggered = triggered & \ - ~np.array(fin['multiple_triggers'][:, trigger_names_dict[indiv_triggers]], dtype=np.bool) + ~np.array(fin['multiple_triggers'][:, trigger_names_dict[indiv_triggers]], dtype=bool) else: for indiv_trigger in indiv_triggers: triggered = triggered & \ - ~np.array(fin['multiple_triggers'][:, trigger_names_dict[indiv_trigger]], dtype=np.bool) + ~np.array(fin['multiple_triggers'][:, trigger_names_dict[indiv_trigger]], dtype=bool) if 'min_sigma' in values.keys(): if isinstance(values['min_sigma'], list): @@ -396,7 +396,7 @@ def get_Veff_Aeff_single( As = np.array(fin['max_amp_ray_solution']) max_amps = np.argmax(As[:, values['ray_channel']], axis=-1) sol = np.array(fin['ray_tracing_solution_type']) - mask = np.array([sol[i, values['ray_channel'], max_amps[i]] == values['ray_solution'] for i in range(len(max_amps))], dtype=np.bool) + mask = np.array([sol[i, values['ray_channel'], max_amps[i]] == values['ray_solution'] for i in range(len(max_amps))], dtype=bool) triggered = triggered & mask if 'n_reflections' in values.keys(): diff --git a/NuRadioMC/utilities/merge_hdf5.py b/NuRadioMC/utilities/merge_hdf5.py index cb31765d9..464f50733 100644 --- a/NuRadioMC/utilities/merge_hdf5.py +++ b/NuRadioMC/utilities/merge_hdf5.py @@ -266,7 +266,7 @@ def merge2(filenames, output_filename): # try: input_files = np.array(sorted(glob.glob(filename + '.part????'))) input_files = np.append(input_files, np.array(sorted(glob.glob(filename + '.part??????')))) - mask = np.array([os.path.getsize(x) > 1000 for x in input_files], dtype=np.bool) + mask = np.array([os.path.getsize(x) > 1000 for x in input_files], dtype=bool) if(np.sum(~mask)): logger.warning("{:d} files were deselected because their filesize was to small".format(np.sum(~mask))) input_args.append({'filenames': input_files[mask], 'output_filename': output_filename}) diff --git a/NuRadioReco/detector/ARA/ara_detector_db.json b/NuRadioReco/detector/ARA/ara_detector_db.json index ec5fa08bb..b0f194bb9 100644 --- a/NuRadioReco/detector/ARA/ara_detector_db.json +++ b/NuRadioReco/detector/ARA/ara_detector_db.json @@ -1,4 +1,3 @@ -{ { "_default": {}, "channels": { diff --git a/NuRadioReco/detector/amp.py b/NuRadioReco/detector/amp.py index 8133af0d6..04dadc714 100644 --- a/NuRadioReco/detector/amp.py +++ b/NuRadioReco/detector/amp.py @@ -25,7 +25,7 @@ def get_amp_response(frequencies, amp_name): get_phase = intp.interp1d(freqs2, np.unwrap(phase)) get_linmag = intp.interp1d(freqs, linmag) - amp_response = np.zeros_like(frequencies, dtype=np.complex) + amp_response = np.zeros_like(frequencies, dtype=complex) mask = (frequencies > max(freqs.min(), freqs2.min())) & (frequencies < min(freqs.max(), freqs2.max())) amp_response[mask] = get_linmag(frequencies[mask]) * np.exp(1j * get_phase(frequencies[mask])) return amp_response diff --git a/NuRadioReco/detector/antennapattern.py b/NuRadioReco/detector/antennapattern.py index 89a6d9795..2452428e1 100644 --- a/NuRadioReco/detector/antennapattern.py +++ b/NuRadioReco/detector/antennapattern.py @@ -64,7 +64,7 @@ def interpolate_linear_vectorized(x, x0, x1, y0, y1, interpolation_method='compl """ x = np.array(x) mask = x0 != x1 - result = np.zeros_like(x, dtype=np.complex) + result = np.zeros_like(x, dtype=complex) denominator = x1 - x0 if interpolation_method == 'complex': result[mask] = y0[mask] + (y1[mask] - y0[mask]) * (x[mask] - x0[mask]) / denominator[mask] @@ -623,13 +623,13 @@ def parse_AERA_XML_file(path): # get frequencies and angles frequencies_node = root.find("./frequency") - frequencies = np.array(frequencies_node.text.strip().split(), dtype=np.float) * units.MHz + frequencies = np.array(frequencies_node.text.strip().split(), dtype=float) * units.MHz theta_node = root.find("./theta") - thetas = np.array(theta_node.text.strip().split(), dtype=np.float) * units.deg + thetas = np.array(theta_node.text.strip().split(), dtype=float) * units.deg phi_node = root.find("./phi") - phis = np.array(phi_node.text.strip().split(), dtype=np.float) * units.deg + phis = np.array(phi_node.text.strip().split(), dtype=float) * units.deg n_freqs = len(frequencies) n_angles = len(phis) @@ -650,16 +650,16 @@ def parse_AERA_XML_file(path): freq_string = "%.1f" % freq theta_amp_node = root.find("./EAHTheta_amp[@idfreq='%s']" % freq_string) - theta_amps[iFreq] = np.array(theta_amp_node.text.strip().split(), dtype=np.float) * units.m + theta_amps[iFreq] = np.array(theta_amp_node.text.strip().split(), dtype=float) * units.m theta_phase_node = root.find("./EAHTheta_phase[@idfreq='%s']" % freq_string) - theta_phases[iFreq] = np.deg2rad(np.array(theta_phase_node.text.strip().split(" "), dtype=np.float)) + theta_phases[iFreq] = np.deg2rad(np.array(theta_phase_node.text.strip().split(" "), dtype=float)) phi_amp_node = root.find("./EAHPhi_amp[@idfreq='%s']" % freq_string) - phi_amps[iFreq] = np.array(phi_amp_node.text.strip().split(), dtype=np.float) * units.m + phi_amps[iFreq] = np.array(phi_amp_node.text.strip().split(), dtype=float) * units.m phi_phase_node = root.find("./EAHPhi_phase[@idfreq='%s']" % freq_string) - phi_phases[iFreq] = np.deg2rad(np.array(phi_phase_node.text.strip().split(), dtype=np.float)) + phi_phases[iFreq] = np.deg2rad(np.array(phi_phase_node.text.strip().split(), dtype=float)) return frequencies, phis, thetas, phi_amps, phi_phases, theta_amps, theta_phases @@ -1059,8 +1059,8 @@ def get_antenna_response_vectorized(self, freq, zenith, azimuth, orientation_the of the same length as the frequency input """ if self._notfound: - VEL = {'theta': np.ones(len(freq), dtype=np.complex), - 'phi': np.ones(len(freq), dtype=np.complex)} + VEL = {'theta': np.ones(len(freq), dtype=complex), + 'phi': np.ones(len(freq), dtype=complex)} return VEL if isinstance(freq, (float, int)): diff --git a/NuRadioReco/detector/filterresponse.py b/NuRadioReco/detector/filterresponse.py index 66dbd1d79..58008e251 100644 --- a/NuRadioReco/detector/filterresponse.py +++ b/NuRadioReco/detector/filterresponse.py @@ -24,7 +24,7 @@ def get_filter_response_mini_circuits(frequencies, filter_name): get_S21 = intp.interp1d(ff, S21) - response = np.zeros_like(frequencies, dtype=np.complex) + response = np.zeros_like(frequencies, dtype=complex) mask = (frequencies > ff.min()) & (frequencies < ff.max()) response[mask] = get_S21(frequencies[mask]) return response @@ -54,7 +54,7 @@ def get_filter_response_mini_circuits2(frequencies, filter_name): phase2 = -2 * np.pi * np.cumsum(get_group_delay(fff2) * df) get_phase = intp.interp1d(fff2, phase2) - response = np.zeros_like(frequencies, dtype=np.complex) + response = np.zeros_like(frequencies, dtype=complex) mask = (frequencies > max(ff.min(), ff2.min())) & (frequencies < min(ff.max(), ff2.max())) response[mask] = get_insertion_loss(frequencies[mask]) * np.exp(1j * get_phase(frequencies[mask])) return response @@ -83,7 +83,7 @@ def get_filter_response(frequencies, filter_name): get_phase = intp.interp1d(ff2, np.unwrap(phase)) get_insertion_loss = intp.interp1d(ff, insertion_loss) - response = np.zeros_like(frequencies, dtype=np.complex) + response = np.zeros_like(frequencies, dtype=complex) mask = (frequencies > max(ff.min(), ff2.min())) & (frequencies < min(ff.max(), ff2.max())) response[mask] = get_insertion_loss(frequencies[mask]) * np.exp(1j * get_phase(frequencies[mask])) return response diff --git a/NuRadioReco/examples/AliasPhasedArray/SNR_study/T01_generate_events_simple.py b/NuRadioReco/examples/AliasPhasedArray/SNR_study/T01_generate_events_simple.py index 9d0947086..8cc768dd2 100644 --- a/NuRadioReco/examples/AliasPhasedArray/SNR_study/T01_generate_events_simple.py +++ b/NuRadioReco/examples/AliasPhasedArray/SNR_study/T01_generate_events_simple.py @@ -47,7 +47,7 @@ vertex_angle_max = 145 * units.deg vertex_angles = np.random.uniform(vertex_angle_min, vertex_angle_max, n_events) -data_sets["yy"] = np.zeros(n_events, dtype=np.float) +data_sets["yy"] = np.zeros(n_events, dtype=float) data_sets["zz"] = z_ant + distance * np.cos(vertex_angles) data_sets["xx"] = distance * np.sin(vertex_angles) @@ -59,7 +59,7 @@ data_sets["event_group_ids"] = np.arange(n_events) data_sets["shower_ids"] = np.arange(n_events) data_sets["n_interaction"] = np.ones(n_events, dtype=int) -data_sets["vertex_times"] = np.zeros(n_events, dtype=np.float) +data_sets["vertex_times"] = np.zeros(n_events, dtype=float) data_sets["flavors"] = np.ones(n_events, dtype=int) * flavor data_sets["energies"] = np.ones(n_events, dtype=int) * energy data_sets["interaction_type"] = ['cc'] * n_events diff --git a/NuRadioReco/examples/RNO_data/read_data_example/read_rnog.py b/NuRadioReco/examples/RNO_data/read_data_example/read_rnog.py new file mode 100644 index 000000000..2855cbe7e --- /dev/null +++ b/NuRadioReco/examples/RNO_data/read_data_example/read_rnog.py @@ -0,0 +1,62 @@ +from NuRadioReco.modules.io.RNO_G import readRNOGDataMattak +from NuRadioReco.modules.io import eventWriter +from NuRadioReco.utilities import units + +import sys +import logging + +""" read in data """ +list_of_root_files = sys.argv[1:-1] +output_filename = sys.argv[-1] + +rnog_reader = readRNOGDataMattak.readRNOGData(log_level=logging.DEBUG) +writer = eventWriter.eventWriter() + +""" +With a selector you can select or reject events based on information in the +Mattak class EventInfo. See https://github.com/RNO-G/mattak/blob/main/py/mattak/Dataset.py + +class EventInfo: + eventNumber: int + station : int + run: int + readoutTime : float + triggerTime : float + triggerType: str + sysclk: int + sysclkLastPPS: Tuple[int, int] # the last 2 PPS sysclks, most recent first + pps: int + radiantStartWindows: numpy.ndarray + sampleRate: float # Sample rate, in GSa/s +""" + +# The following selector selects only events with a forced trigger. +selectors = [lambda einfo: einfo.triggerType == "FORCE"] + +rnog_reader.begin( + list_of_root_files, + selectors=selectors, + # Currently false because Mattak does not contain calibrated data yet + read_calibrated_data=False, + # Only used when read_calibrated_data==False, performs a simple baseline subtraction each 128 bins + apply_baseline_correction=True, + # Only used when read_calibrated_data==False, performs a linear voltage calibration with hardcoded values + convert_to_voltage=True, + # Can be used instead of defining a selector (only for triggers) + select_triggers=None, + # If true, and if the RunTable database is available select runs based on the following criteria + select_runs=True, + # Only use runs of a certain run type + run_types=["physics"], + # Only use runs with a maximum trigger rate of 1 Hz + max_trigger_rate=1 * units.Hz) + +writer.begin(filename=output_filename) + +for i_event, event in enumerate(rnog_reader.run()): + writer.run(event) + +rnog_reader.end() +writer.end() + + diff --git a/NuRadioReco/examples/cr_efficiency_analysis/helper_cr_eff.py b/NuRadioReco/examples/cr_efficiency_analysis/helper_cr_eff.py index 244f1a916..0545f2954 100644 --- a/NuRadioReco/examples/cr_efficiency_analysis/helper_cr_eff.py +++ b/NuRadioReco/examples/cr_efficiency_analysis/helper_cr_eff.py @@ -13,9 +13,9 @@ class NumpyEncoder(json.JSONEncoder): """ Special json encoder for numpy types """ def default(self, obj): - if isinstance(obj, np.integer): + if isinstance(obj, int): return int(obj) - elif isinstance(obj, np.floating): + elif isinstance(obj, float): return float(obj) elif isinstance(obj, np.ndarray): return obj.tolist() diff --git a/NuRadioReco/framework/base_trace.py b/NuRadioReco/framework/base_trace.py index 0b735967f..c79c86c29 100644 --- a/NuRadioReco/framework/base_trace.py +++ b/NuRadioReco/framework/base_trace.py @@ -36,7 +36,7 @@ def get_trace(self): trace: np.array of floats the time trace """ - if(not self.__time_domain_up_to_date): + if not self.__time_domain_up_to_date: self._time_trace = fft.freq2time(self._frequency_spectrum, self._sampling_rate) self.__time_domain_up_to_date = True self._frequency_spectrum = None @@ -62,10 +62,10 @@ def get_filtered_trace(self, passband, filter_type='butter', order=10, rp=None): return fft.freq2time(spec, self.get_sampling_rate()) def get_frequency_spectrum(self): - if(self.__time_domain_up_to_date): + if self.__time_domain_up_to_date: self._frequency_spectrum = fft.time2freq(self._time_trace, self._sampling_rate) self._time_trace = None -# logger.debug("frequency spectrum has shape {}".format(self._frequency_spectrum.shape)) + # logger.debug("frequency spectrum has shape {}".format(self._frequency_spectrum.shape)) self.__time_domain_up_to_date = False return np.copy(self._frequency_spectrum) @@ -82,7 +82,8 @@ def set_trace(self, trace, sampling_rate): """ if trace is not None: if trace.shape[trace.ndim - 1] % 2 != 0: - raise ValueError('Attempted to set trace with an uneven number ({}) of samples. Only traces with an even number of samples are allowed.'.format(trace.shape[trace.ndim - 1])) + raise ValueError(('Attempted to set trace with an uneven number ({}) of samples. ' + 'Only traces with an even number of samples are allowed.').format(trace.shape[trace.ndim - 1])) self.__time_domain_up_to_date = True self._time_trace = np.copy(trace) self._sampling_rate = sampling_rate @@ -109,9 +110,10 @@ def get_times(self): try: length = self.get_number_of_samples() times = np.arange(0, length / self._sampling_rate - 0.1 / self._sampling_rate, 1. / self._sampling_rate) + self._trace_start_time - if(len(times) != length): - logger.error("time array does not have the same length as the trace. n_samples = {:d}, sampling rate = {:.5g}".format(length, self._sampling_rate)) - raise ValueError("time array does not have the same length as the trace") + if len(times) != length: + err = f"time array does not have the same length as the trace. n_samples = {length:d}, sampling rate = {self._sampling_rate:.5g}" + logger.error(err) + raise ValueError(err) except: times = np.array([]) return times @@ -148,7 +150,7 @@ def get_number_of_samples(self): n_samples: int number of samples in time domain """ - if(self.__time_domain_up_to_date): + if self.__time_domain_up_to_date: length = self._time_trace.shape[-1] # returns the correct length independent of the dimension of the array (channels are 1dim, efields are 3dim) else: length = (self._frequency_spectrum.shape[-1] - 1) * 2 @@ -184,12 +186,14 @@ def resample(self, sampling_rate): if resampling_factor.numerator != 1: # resample and use axis -1 since trace might be either shape (N) for analytic trace or shape (3,N) for E-field resampled_trace = scipy.signal.resample(resampled_trace, resampling_factor.numerator * self.get_number_of_samples(), axis=-1) + if resampling_factor.denominator != 1: # resample and use axis -1 since trace might be either shape (N) for analytic trace or shape (3,N) for E-field resampled_trace = scipy.signal.resample(resampled_trace, np.shape(resampled_trace)[-1] // resampling_factor.denominator, axis=-1) if resampled_trace.shape[-1] % 2 != 0: resampled_trace = resampled_trace.T[:-1].T + self.set_trace(resampled_trace, sampling_rate) def serialize(self): @@ -201,7 +205,7 @@ def serialize(self): def deserialize(self, data_pkl): data = pickle.loads(data_pkl) self.set_trace(data['time_trace'], data['sampling_rate']) - if('trace_start_time' in data.keys()): + if 'trace_start_time' in data.keys(): self.set_trace_start_time(data['trace_start_time']) def __add__(self, x): @@ -214,10 +218,13 @@ def __add__(self, x): # Some sanity checks if not isinstance(x, BaseTrace): raise TypeError('+ operator is only defined for 2 BaseTrace objects') + if self.get_trace() is None or x.get_trace() is None: raise ValueError('One of the trace objects has no trace set') + if self.get_trace().ndim != x.get_trace().ndim: raise ValueError('Traces have different dimensions') + if self.get_sampling_rate() != x.get_sampling_rate(): # Upsample trace with lower sampling rate # Create new baseTrace object for the resampling so we don't change the originals @@ -249,10 +256,12 @@ def __add__(self, x): first_trace = trace_2 second_trace = trace_1 trace_start = x.get_trace_start_time() + # Calculate the difference in the trace start time between the traces and the number of # samples that time difference corresponds to time_offset = np.abs(x.get_trace_start_time() - self.get_trace_start_time()) i_start = int(round(time_offset * sampling_rate)) + # We have to distinguish 2 cases: Trace is 1D (channel) or 2D(E-field) # and treat them differently if trace_1.ndim == 1: @@ -273,11 +282,13 @@ def __add__(self, x): early_trace[:, :first_trace.shape[1]] = first_trace late_trace = np.zeros((second_trace.shape[0], trace_length)) late_trace[:, :second_trace.shape[1]] = second_trace + # Correct for different trace start times by using fourier shift theorem to # shift the later trace backwards. late_trace_object = BaseTrace() late_trace_object.set_trace(late_trace, sampling_rate) late_trace_object.apply_time_shift(time_offset, True) + # Create new BaseTrace object holding the summed traces new_trace = BaseTrace() new_trace.set_trace(early_trace + late_trace_object.get_trace(), sampling_rate) diff --git a/NuRadioReco/modules/analogToDigitalConverter.py b/NuRadioReco/modules/analogToDigitalConverter.py index 37b229054..528666cf8 100644 --- a/NuRadioReco/modules/analogToDigitalConverter.py +++ b/NuRadioReco/modules/analogToDigitalConverter.py @@ -50,7 +50,7 @@ def perfect_comparator(trace, adc_n_bits, adc_ref_voltage, mode='floor', output= digital_trace = round_to_int(digital_trace) if (output == 'voltage'): - digital_trace = lsb_voltage * digital_trace.astype(np.float) + digital_trace = lsb_voltage * digital_trace.astype(float) elif (output == 'counts'): pass else: diff --git a/NuRadioReco/modules/channelGalacticNoiseAdder.py b/NuRadioReco/modules/channelGalacticNoiseAdder.py index 0638fb423..9afafc989 100644 --- a/NuRadioReco/modules/channelGalacticNoiseAdder.py +++ b/NuRadioReco/modules/channelGalacticNoiseAdder.py @@ -122,7 +122,7 @@ def run( passband_filter = (freqs > passband[0]) & (freqs < passband[1]) noise_spec_sum = np.zeros_like(channel.get_frequency_spectrum()) flux_sum = np.zeros(freqs[passband_filter].shape) - efield_sum = np.zeros((3, freqs.shape[0]), dtype=np.complex) + efield_sum = np.zeros((3, freqs.shape[0]), dtype=complex) if self.__debug: plt.close('all') fig = plt.figure(figsize=(12, 8)) @@ -177,7 +177,7 @@ def run( ax_3.plot(freqs[passband_filter] / units.MHz, E / (units.V / units.m), c='k', alpha=.02) # assign random phases and polarizations to electric field - noise_spectrum = np.zeros((3, freqs.shape[0]), dtype=np.complex) + noise_spectrum = np.zeros((3, freqs.shape[0]), dtype=complex) phases = np.random.uniform(0, 2. * np.pi, len(S)) polarizations = np.random.uniform(0, 2. * np.pi, len(S)) diff --git a/NuRadioReco/modules/channelGenericNoiseAdder.py b/NuRadioReco/modules/channelGenericNoiseAdder.py index e8a4acfd9..e5556c02a 100644 --- a/NuRadioReco/modules/channelGenericNoiseAdder.py +++ b/NuRadioReco/modules/channelGenericNoiseAdder.py @@ -2,7 +2,7 @@ from NuRadioReco.modules.base.module import register_run import numpy as np from NuRadioReco.utilities import units, fft -import numpy.random +from numpy.random import Generator, Philox import logging @@ -27,7 +27,7 @@ def add_random_phases(self, amps, n_samples_time_domain): """ amps = np.array(amps, dtype='complex') Np = (n_samples_time_domain - 1) // 2 - phases = self.__random_generator.rand(Np) * 2 * np.pi + phases = self.__random_generator.random(Np) * 2 * np.pi phases = np.cos(phases) + 1j * np.sin(phases) amps[1:Np + 1] *= phases # Note that the last entry of the index slice is f[Np] ! @@ -45,7 +45,7 @@ def fftnoise_fullfft(self, f): """ f = np.array(f, dtype='complex') Np = (len(f) - 1) // 2 - phases = self.__random_generator.rand(Np) * 2 * np.pi + phases = self.__random_generator.random(Np) * 2 * np.pi phases = np.cos(phases) + 1j * np.sin(phases) f[1:Np + 1] *= phases # Note that the last entry of the index slice is f[Np] ! f[-1:-1 - Np:-1] = np.conj(f[1:Np + 1]) @@ -308,7 +308,7 @@ def __init__(self): def begin(self, debug=False, seed=None): self.__debug = debug - self.__random_generator = np.random.RandomState(seed) + self.__random_generator = Generator(Philox(seed)) if debug: self.logger.setLevel(logging.DEBUG) diff --git a/NuRadioReco/modules/electricFieldSignalReconstructor.py b/NuRadioReco/modules/electricFieldSignalReconstructor.py index 94252131b..9f4089d2b 100644 --- a/NuRadioReco/modules/electricFieldSignalReconstructor.py +++ b/NuRadioReco/modules/electricFieldSignalReconstructor.py @@ -98,9 +98,9 @@ def run(self, evt, station, det, debug=False): times = electric_field.get_times() 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=np.bool) + 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=np.bool) # the last n bins + 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 signal_energy_fluence = trace_utilities.get_electric_field_energy_fluence(trace, times, mask_signal_window, mask_noise_window) dt = times[1] - times[0] diff --git a/NuRadioReco/modules/envelope_phasedarray/triggerSimulator.py b/NuRadioReco/modules/envelope_phasedarray/triggerSimulator.py index 0a7aad5fa..e3df099c8 100644 --- a/NuRadioReco/modules/envelope_phasedarray/triggerSimulator.py +++ b/NuRadioReco/modules/envelope_phasedarray/triggerSimulator.py @@ -102,7 +102,7 @@ def envelope_trigger(self, adc_type='perfect_floor_comparator', diode=diode) time_step = 1 / det.get_channel(station_id, channel_id)['trigger_adc_sampling_frequency'] - times = np.arange(len(trace), dtype=np.float) * time_step + times = np.arange(len(trace), dtype=float) * time_step times += channel.get_trace_start_time() else: diff --git a/NuRadioReco/modules/io/rno_g/__init__.py b/NuRadioReco/modules/io/RNO_G/__init__.py similarity index 100% rename from NuRadioReco/modules/io/rno_g/__init__.py rename to NuRadioReco/modules/io/RNO_G/__init__.py diff --git a/NuRadioReco/modules/io/RNO_G/readRNOGDataMattak.py b/NuRadioReco/modules/io/RNO_G/readRNOGDataMattak.py new file mode 100644 index 000000000..db2835ecb --- /dev/null +++ b/NuRadioReco/modules/io/RNO_G/readRNOGDataMattak.py @@ -0,0 +1,792 @@ +import numpy as np +import logging +import os +import time +import astropy.time +import math + +from NuRadioReco.modules.base.module import register_run + +import NuRadioReco.framework.event +import NuRadioReco.framework.station +import NuRadioReco.framework.channel +import NuRadioReco.framework.trigger + +from NuRadioReco.utilities import units +import mattak.Dataset + + +def baseline_correction(wfs, n_bins=128, func=np.median): + """ + Simple baseline correction function. Determines baseline in discrete chuncks of "n_bins" with + the function specified (i.e., mean or median). + + Parameters + ---------- + + wfs: np.array(n_events, n_channels, n_samples) + Waveforms of several events/channels. + + n_bins: int + Number of samples/bins in one "chunck". If None, calculate median/mean over entire trace. (Default: 128) + + func: np.mean or np.median + Function to calculate pedestal + + Returns + ------- + + wfs_corrected: np.array(n_events, n_channels, n_samples) + Baseline/pedestal corrected waveforms + """ + + # Example: Get baselines in chunks of 128 bins + # wfs in (n_events, n_channels, 2048) + # np.split -> (16, n_events, n_channels, 128) each waveform split in 16 chuncks + # func -> (16, n_events, n_channels) pedestal for each chunck + if n_bins is not None: + baseline_values = func(np.split(wfs, 2048 // n_bins, axis=-1), axis=-1) + + # np.repeat -> (2048, n_events, n_channels) concatenate the 16 chuncks to one baseline + baseline_traces = np.repeat(baseline_values, n_bins % 2048, axis=0) + else: + baseline_values = func(wfs, axis=-1) + # np.repeat -> (2048, n_events, n_channels) concatenate the 16 chuncks to one baseline + baseline_traces = np.repeat(baseline_values, 2048, axis=0) + + # np.moveaxis -> (n_events, n_channels, 2048) + baseline_traces = np.moveaxis(baseline_traces, 0, -1) + + return wfs - baseline_traces + + +def get_time_offset(trigger_type): + """ + Mapping the offset between trace start time and trigger time (~ signal time). + Temporary use hard-coded values for each trigger type. In the future this + information might be time, station, and channel dependent and should come + from a database (or is already calibrated in mattak) + + Current values motivated by figures posted in PR https://github.com/nu-radio/NuRadioMC/pull/519 + + Parameters + ---------- + + trigger_type: str + Trigger type encoded as string from Mattak + + Returns + ------- + + time_offset: float + trace_start_time = trigger_time - time_offset + + """ + + time_offsets = { + "FORCE": 0, + "LT": 250 * units.ns, + "RADIANT": 475 * units.ns, + "UNKNOWN": 0 # Due to a firmware issue at the beginning of data taking the trigger types were not properly set. + } + + # Should have the same time offset ?! + if trigger_type.startswith("RADIANT"): + trigger_type = "RADIANT" + + if trigger_type in time_offsets: + return time_offsets[trigger_type] + else: + known_trigger_types = ", ".join(time_offsets.keys()) + raise KeyError(f"Unknown trigger type: {trigger_type}. Known are: {known_trigger_types}. Abort ....") + + +def all_files_in_directory(mattak_dir): + """ + Checks if all Mattak root files are in a directory. + Ignoring runinfo.root because (asaik) not all runs have those and information is currently not read by Mattak. + There are mattak directories which produce a ReferenceError when reading. They have a "combined.root" which is + apparently empty but are missing the daqstatus, pedestal, and header file. + + Parameters + ---------- + + mattak_dir: str + Path to a mattak directory + + Returns + ------- + + all_there: bool + True, if all "req_files" are there and waveforms.root or combined.root. Otherwise returns False. + """ + # one or the other has to be present + if not os.path.exists(os.path.join(mattak_dir, "waveforms.root")) and \ + not os.path.exists(os.path.join(mattak_dir, "combined.root")): + return False + + req_files = ["daqstatus.root", "headers.root", "pedestal.root"] + for file in req_files: + if not os.path.exists(os.path.join(mattak_dir, file)): + return False + + return True + + +class readRNOGData: + + def __init__(self, run_table_path=None, log_level=logging.INFO): + """ + Parameters + ---------- + + run_table_path: str + Path to a run_table.cvs file. If None, the run table is queried from the DB. (Default: None) + + log_level: enum + Set verbosity level of logger. If logging.DEBUG, set mattak to verbose (unless specified in mattak_kwargs). + (Default: logging.INFO) + """ + self.logger = logging.getLogger('NuRadioReco.readRNOGData') + self.logger.setLevel(log_level) + + # Initialize run table for run selection + self.__run_table = None + + if run_table_path is None: + try: + from rnog_data.runtable import RunTable + self.logger.debug("Access RunTable database ...") + try: + self.__run_table = RunTable().get_table() + except: + self.logger.warn("No connect to RunTable database could be established. " + "Runs can not be filtered.") + except ImportError: + self.logger.warn("Import of run table failed. Runs can not be filtered.! \n" + "You can get the interface from GitHub: git@github.com:RNO-G/rnog-runtable.git") + else: + import pandas + self.__run_table = pandas.read_csv(run_table_path) + + + def begin(self, + dirs_files, + read_calibrated_data=False, + select_triggers=None, + select_runs=False, + apply_baseline_correction=True, + convert_to_voltage=True, + selectors=[], + run_types=["physics"], + run_time_range=None, + max_trigger_rate=0 * units.Hz, + mattak_kwargs={}, + overwrite_sampling_rate=None): + """ + Parameters + ---------- + + dirs_files: str, list(str) + Path to run directories (i.e. ".../stationXX/runXXX/") or path to root files (have to be "combined" mattak files). + + read_calibrated_data: bool + If True, read calibrated waveforms from Mattak.Dataset. If False, read "raw" ADC traces. + (temp. Default: False, this can/should be switched once the calibration in incorp. into Mattak) + + select_triggers: str or list(str) + Names of triggers which should be selected. Convinence interface instead of passing a selector + (see "selectors" below). (Default: None) + + select_runs: bool + If True, use information in run_table to select runs (based on run_type, run_time, trigger_rate, ...). + If the run_table is not available no selection is performed (and the programm is not interrupted, + only an error message is raised). See parameters to configure run selection. (Default: False) + + Other Parameters + ---------------- + + apply_baseline_correction: bool + Only applies when non-calibrated data are read. If true, correct for DC offset. + (Default: True) + + convert_to_voltage: bool + Only applies when non-calibrated data are read. If true, convert ADC to voltage. + (Default: True) + + selectors: list of lambdas + List of lambda(eventInfo) -> bool to pass to mattak.Dataset.iterate to select events. + Example: trigger_selector = lambda eventInfo: eventInfo.triggerType == "FORCE" + + run_types: list + Used to select/reject runs from information in the RNO-G RunTable. List of run_types to be used. (Default: ['physics']) + + run_time_range: tuple + Specify a time range to select runs (it is sufficient that runs cover the time range partially). + Each value of the tuple has to be in a format which astropy.time.Time understands. A value can be None + which means that the lower or upper bound is unconstrained. If run_time_range is None no time selection is + applied. (Default: None) + + max_trigger_rate: float + Used to select/reject runs from information in the RNO-G RunTable. Maximum allowed trigger rate (per run) in Hz. + If 0, no cut is applied. (Default: 1 Hz) + + mattak_kwargs: dict + Dictionary of arguments for mattak.Dataset.Dataset. (Default: {}) + Example: Select a mattak "backend". Options are "auto", "pyroot", "uproot". If "auto" is selected, + pyroot is used if available otherwise a "fallback" to uproot is used. (Default: "auto") + + overwrite_sampling_rate: float + Set sampling rate of the imported waveforms. This overwrites what is read out from runinfo (i.e., stored in the mattak files). + If None, nothing is overwritten and the sampling rate from the mattak file is used. (Default: None) + NOTE: This option might be necessary when old mattak files are read which have this not set. + """ + + t0 = time.time() + + self._read_calibrated_data = read_calibrated_data + self._apply_baseline_correction = apply_baseline_correction + self._convert_to_voltage = convert_to_voltage + + # Temporary solution hard-coded values from Cosmin. Only used when uncalibrated data + # is read and convert_to_voltage is True. + self._adc_ref_voltage_range = 2.5 * units.volt + self._adc_n_bits = 12 + + self._overwrite_sampling_rate = overwrite_sampling_rate + + # Set parameter for run selection + self.__max_trigger_rate = max_trigger_rate + self.__run_types = run_types + + if run_time_range is not None: + convert_time = lambda t: None if t is None else astropy.time.Time(t) + self._time_low = convert_time(run_time_range[0]) + self._time_high = convert_time(run_time_range[1]) + else: + self._time_low = None + self._time_high = None + + if select_runs and self.__run_table is not None: + self.logger.info("\n\tSelect runs with type: {}".format(", ".join(run_types)) + + f"\n\tSelect runs with max. trigger rate of {max_trigger_rate / units.Hz} Hz" + f"\n\tSelect runs which are between {self._time_low} - {self._time_high}") + + self.set_selectors(selectors, select_triggers) + + # Read data + self._time_begin = 0 + self._time_run = 0 + self.__counter = 0 + self.__skipped = 0 + self.__invalid = 0 + + self._events_information = None + self._datasets = [] + self.__n_events_per_dataset = [] + + self.logger.info(f"Parse through / read-in {len(dirs_files)} directory(ies) / file(s).") + + self.__skipped_runs = 0 + self.__n_runs = 0 + + if not isinstance(dirs_files, (list, np.ndarray)): + dirs_files = [dirs_files] + + # Set verbose for mattak + if "verbose" in mattak_kwargs: + verbose = mattak_kwargs.pop("verbose") + else: + verbose = self.logger.level >= logging.DEBUG + + for dir_file in dirs_files: + + if not os.path.exists(dir_file): + self.logger.error(f"The directory/file {dir_file} does not exist") + continue + + if os.path.isdir(dir_file): + + if not all_files_in_directory(dir_file): + self.logger.error(f"Incomplete directory: {dir_file}. Skip ...") + continue + + try: + dataset = mattak.Dataset.Dataset(station=0, run=0, data_dir=dir_file, verbose=verbose, **mattak_kwargs) + except (ReferenceError, KeyError) as e: + self.logger.error(f"The following exeption was raised reading in the run: {dir_file}. Skip that run ...:\n", exc_info=e) + continue + else: + raise NotImplementedError("The option to read in files is not implemented yet") + + + # filter runs/datasets based on + if select_runs and self.__run_table is not None and not self.__select_run(dataset): + self.__skipped_runs += 1 + continue + + self.__n_runs += 1 + self._datasets.append(dataset) + self.__n_events_per_dataset.append(dataset.N()) + + if not len(self._datasets): + err = "Found no valid datasets. Stop!" + self.logger.error(err) + raise FileNotFoundError(err) + + # keeps track which event index is in which dataset + self._event_idxs_datasets = np.cumsum(self.__n_events_per_dataset) + self._n_events_total = np.sum(self.__n_events_per_dataset) + self._time_begin = time.time() - t0 + + self.logger.info(f"{self._n_events_total} events in {len(self._datasets)} runs/datasets " + f"have been found using the {self._datasets[0].backend} Mattak backend.") + + if not self._n_events_total: + err = "No runs have been selected. Abort ..." + self.logger.error(err) + raise ValueError(err) + + + def set_selectors(self, selectors, select_triggers=None): + """ + Parameters + ---------- + + selectors: list of lambdas + List of lambda(eventInfo) -> bool to pass to mattak.Dataset.iterate to select events. + Example: trigger_selector = lambda eventInfo: eventInfo.triggerType == "FORCE" + + select_triggers: str or list(str) + Names of triggers which should be selected. Convinence interface instead of passing a selector. (Default: None) + """ + + # Initialize selectors for event filtering + if not isinstance(selectors, (list, np.ndarray)): + selectors = [selectors] + + if select_triggers is not None: + if isinstance(select_triggers, str): + selectors.append(lambda eventInfo: eventInfo.triggerType == select_triggers) + else: + for select_trigger in select_triggers: + selectors.append(lambda eventInfo: eventInfo.triggerType == select_trigger) + + self._selectors = selectors + self.logger.info(f"Set {len(self._selectors)} selector(s)") + + + def __select_run(self, dataset): + """ Filter/select runs/datasets. + + Parameters + ---------- + + dataset: mattak.Dataset.Dataset + + select: bool + Return True to select an dataset, return False to reject/skip it. + """ + + # get first eventInfo + dataset.setEntries(0) + event_info = dataset.eventInfo() + + run_id = event_info.run + station_id = event_info.station + + run_info = self.__run_table.query(f"station == {station_id:d} & run == {run_id:d}") + + if not len(run_info): + self.logger.error(f"Run {run_id:d} (station {station_id:d}) not in run table. Reject...") + return False + + # "time_start/end" is stored in the isot format. datetime is much faster than astropy (~85ns vs 55 mus). + # But using datetime would mean to stip decimals because datetime can only handle mu sec precision and can not cope + # with the additional decimals for ns. + if self._time_low is not None: + time_end = astropy.time.Time(run_info["time_end"].values[0]) + if time_end < self._time_low: + self.logger.info(f"Reject station {station_id} run {run_id} because run ended before {self._time_low}") + return False + + if self._time_high is not None: + time_start = astropy.time.Time(run_info["time_start"].values[0]) + if time_start > self._time_high: + self.logger.info(f"Reject station {station_id} run {run_id} because run started time after {self._time_high}") + return False + + run_type = run_info["run_type"].values[0] + if not run_type in self.__run_types: + self.logger.info(f"Reject station {station_id} run {run_id} because of run type {run_type}") + return False + + trigger_rate = run_info["trigger_rate"].values[0] * units.Hz + if self.__max_trigger_rate and trigger_rate > self.__max_trigger_rate: + self.logger.info(f"Reject station {station_id} run {run_id} because trigger rate is to high ({trigger_rate / units.Hz:.2f} Hz)") + return False + + return True + + + def __get_n_events_of_prev_datasets(self, dataset_idx): + """ Get accumulated number of events from previous datasets """ + dataset_idx_prev = dataset_idx - 1 + return int(self._event_idxs_datasets[dataset_idx_prev]) if dataset_idx_prev >= 0 else 0 + + + def __get_dataset_for_event(self, event_idx): + """ Get correct dataset and set entry accordingly to event index + + Parameters + ---------- + + event_index: int + Same as in read_event(). + + Returns + ------- + + dataset: mattak.Dataset.Dataset + """ + # find correct dataset + dataset_idx = np.digitize(event_idx, self._event_idxs_datasets) + dataset = self._datasets[dataset_idx] + + event_idx_in_dataset = event_idx - self.__get_n_events_of_prev_datasets(dataset_idx) + dataset.setEntries(event_idx_in_dataset) # increment iterator -> point to new event + + return dataset + + + def _filter_event(self, evtinfo, event_idx=None): + """ Filter an event base on its EventInfo and the configured selectors. + + Parameters + ---------- + + event_info: mattak.Dataset.EventInfo + The event info object for one event. + + event_index: int + Same as in read_event(). Only use for logger.info(). (Default: None) + + Returns + ------- + + skip: bool + Returns True to skip/reject event, return False to keep/read event + """ + if self._selectors is not None: + for selector in self._selectors: + if not selector(evtinfo): + self.logger.debug(f"Event {event_idx} (station {evtinfo.station}, run {evtinfo.run}, " + f"event number {evtinfo.eventNumber}) is skipped.") + self.__skipped += 1 + return True + + return False + + + def get_events_information(self, keys=["station", "run", "eventNumber"]): + """ Return information of all events from the EventInfo + + This function is useful to make a pre-selection of events before actually reading them in combination with + self.read_event(). + + Parameters + ---------- + + keys: list(str) + List of the information to receive from each event. Have to match the attributes (member variables) + of the mattak.Dataset.EventInfo class (examples are "station", "run", "triggerTime", "triggerType", "eventNumber", ...). + (Default: ["station", "run", "eventNumber"]) + + Returns + ------- + + data: dict + Keys of the dict are the event indecies (as used in self.read_event(event_index)). The values are dictinaries + them self containing the information specified with "keys" parameter. + """ + + # Read if dict is None ... + do_read = self._events_information is None + + if not do_read: + # ... or when it does not have the desired information + first_event_info = self._events_information[list(self._events_information.keys())[0]] + + for key in keys: + if key not in list(first_event_info.keys()): + do_read = True + + if do_read: + + self._events_information = {} + n_prev = 0 + for dataset in self._datasets: + dataset.setEntries((0, dataset.N())) + + for idx, evtinfo in enumerate(dataset.eventInfo()): # returns a list + + event_idx = idx + n_prev # event index accross all datasets combined + + if self._filter_event(evtinfo, event_idx): + continue + + self._events_information[event_idx] = {key: getattr(evtinfo, key) for key in keys} + + n_prev += dataset.N() + + return self._events_information + + + def _check_for_valid_information_in_event_info(self, event_info): + """ + Checks if certain information (sampling rate, trigger time) in mattak.Dataset.EventInfo are valid + + Parameters + ---------- + + event_info: mattak.Dataset.EventInfo + + Returns + ------- + + is_valid: bool + Returns True if all information valid, false otherwise + """ + + + if math.isinf(event_info.triggerTime): + self.logger.error(f"Event {event_info.eventNumber} (st {event_info.station}, run {event_info.run}) " + "has inf trigger time. Skip event...") + self.__invalid += 1 + return False + + + if (event_info.sampleRate == 0 or event_info.sampleRate is None) and self._overwrite_sampling_rate is None: + self.logger.error(f"Event {event_info.eventNumber} (st {event_info.station}, run {event_info.run}) " + f"has a sampling rate of {event_info.sampleRate:.2f} GHz. Event is skipped ... " + f"You can avoid this by setting 'overwrite_sampling_rate' in the begin() method.") + self.__invalid += 1 + return False + + return True + + + def _get_event(self, event_info, waveforms): + """ Return a NuRadioReco event + + Parameters + ---------- + + event_info: mattak.Dataset.EventInfo + The event info object for one event. + + waveforms: np.array(n_channel, n_samples) + Typically what dataset.wfs() returns (for one event!) + + Returns + ------- + + evt: NuRadioReco.framework.event + """ + + trigger_time = event_info.triggerTime + if self._overwrite_sampling_rate is not None: + sampling_rate = self._overwrite_sampling_rate + else: + sampling_rate = event_info.sampleRate + + evt = NuRadioReco.framework.event.Event(event_info.run, event_info.eventNumber) + station = NuRadioReco.framework.station.Station(event_info.station) + station.set_station_time(astropy.time.Time(trigger_time, format='unix')) + + trigger = NuRadioReco.framework.trigger.Trigger(event_info.triggerType) + trigger.set_triggered() + trigger.set_trigger_time(trigger_time) + station.set_trigger(trigger) + + for channel_id, wf in enumerate(waveforms): + channel = NuRadioReco.framework.channel.Channel(channel_id) + if self._read_calibrated_data: + channel.set_trace(wf * units.mV, sampling_rate * units.GHz) + else: + # wf stores ADC counts + + if self._apply_baseline_correction: + # correct baseline + wf = baseline_correction(wf) + + if self._convert_to_voltage: + # convert adc to voltage + wf *= (self._adc_ref_voltage_range / (2 ** (self._adc_n_bits) - 1)) + + channel.set_trace(wf, sampling_rate * units.GHz) + + time_offset = get_time_offset(event_info.triggerType) + channel.set_trace_start_time(-time_offset) # relative to event/trigger time + + station.add_channel(channel) + + evt.set_station(station) + + return evt + + + @register_run() + def run(self): + """ + Loop over all events. + + Returns + ------- + + evt: generator(NuRadioReco.framework.event) + """ + event_idx = -1 + for dataset in self._datasets: + dataset.setEntries((0, dataset.N())) + + # read all event infos of the entier dataset (= run) + event_infos = dataset.eventInfo() + wfs = None + + for idx, evtinfo in enumerate(event_infos): # returns a list + event_idx += 1 + + self.logger.debug(f"Processing event number {event_idx} out of total {self._n_events_total}") + t0 = time.time() + + if self._filter_event(evtinfo, event_idx): + continue + + if not self._check_for_valid_information_in_event_info(evtinfo): + continue + + # Just read wfs if necessary + if wfs is None: + wfs = dataset.wfs() + + waveforms_of_event = wfs[idx] + + evt = self._get_event(evtinfo, waveforms_of_event) + + self._time_run += time.time() - t0 + self.__counter += 1 + + yield evt + + + + def get_event_by_index(self, event_index): + """ Allows to read a specific event identifed by its index + + Parameters + ---------- + + event_index: int + The index of a particluar event. The index is the chronological number from 0 to + number of total events (across all datasets). + + Returns + ------- + + evt: NuRadioReco.framework.event + """ + + self.logger.debug(f"Processing event number {event_index} out of total {self._n_events_total}") + t0 = time.time() + + dataset = self.__get_dataset_for_event(event_index) + event_info = dataset.eventInfo() # returns a single eventInfo + + if self._filter_event(event_info, event_index): + return None + + # check this before reading the wfs + if not self._check_for_valid_information_in_event_info(event_info): + return None + + # access data + waveforms = dataset.wfs() + + evt = self._get_event(event_info, waveforms) + + self._time_run += time.time() - t0 + self.__counter += 1 + + return evt + + + def get_event(self, run_nr, event_id): + """ Allows to read a specific event identifed by run number and event id + + Parameters + ---------- + + run_nr: int + Run number + + event_id: int + Event Id + + Returns + ------- + + evt: NuRadioReco.framework.event + """ + + self.logger.debug(f"Processing event {event_id}") + t0 = time.time() + + event_infos = self.get_events_information(keys=["eventNumber", "run"]) + event_idx_ids = np.array([[index, ele["eventNumber"], ele["run"]] for index, ele in event_infos.items()]) + mask = np.all([event_idx_ids[:, 1] == event_id, event_idx_ids[:, 2] == run_nr], axis=0) + + if not np.any(mask): + self.logger.info(f"Could not find event with id: {event_id}.") + return None + elif np.sum(mask) > 1: + self.logger.error(f"Found several events with the same id: {event_id}.") + raise ValueError(f"Found several events with the same id: {event_id}.") + else: + pass + + # int(...) necessary to pass it to mattak + event_index = int(event_idx_ids[mask, 0][0]) + + dataset = self.__get_dataset_for_event(event_index) + event_info = dataset.eventInfo() # returns a single eventInfo + + if self._filter_event(event_info, event_index): + return None + + # check this before reading the wfs + if not self._check_for_valid_information_in_event_info(event_info): + return None + + # access data + waveforms = dataset.wfs() + + evt = self._get_event(event_info, waveforms) + + self._time_run += time.time() - t0 + self.__counter += 1 + + return evt + + + def end(self): + if self.__counter: + self.logger.info( + f"\n\tRead {self.__counter} events (skipped {self.__skipped} events, {self.__invalid} invalid events)" + f"\n\tTime to initialize data sets : {self._time_begin:.2f}s" + f"\n\tTime to read all events : {self._time_run:.2f}s" + f"\n\tTime to per event : {self._time_run / self.__counter:.2f}s" + f"\n\tRead {self.__n_runs} runs, skipped {self.__skipped_runs} runs.") + else: + self.logger.info( + f"\n\tRead {self.__counter} events (skipped {self.__skipped} events, {self.__invalid} invalid events)") \ No newline at end of file diff --git a/NuRadioReco/modules/io/eventWriter.py b/NuRadioReco/modules/io/eventWriter.py index a924b646d..0204943bd 100644 --- a/NuRadioReco/modules/io/eventWriter.py +++ b/NuRadioReco/modules/io/eventWriter.py @@ -72,12 +72,14 @@ def begin(self, filename, max_file_size=1024, check_for_duplicates=False, events both set, the file will be split whenever any of the two conditions is fullfilled. """ logger.setLevel(log_level) - if filename[-4:] == '.nur': + if filename.endswith(".nur"): self.__filename = filename[:-4] else: self.__filename = filename - if filename[-4:] == '.ari': + + if filename.endswith('.ari'): logger.warning('The file ending .ari for NuRadioReco files is deprecated. Please use .nur instead.') + self.__check_for_duplicates = check_for_duplicates self.__number_of_events = 0 self.__current_file_size = 0 diff --git a/NuRadioReco/modules/io/rno_g/readRNOGData.py b/NuRadioReco/modules/io/rno_g/readRNOGData.py deleted file mode 100755 index e1fc23d08..000000000 --- a/NuRadioReco/modules/io/rno_g/readRNOGData.py +++ /dev/null @@ -1,224 +0,0 @@ -import NuRadioReco.framework.event -from NuRadioReco.modules.base.module import register_run -import NuRadioReco.framework.station -import NuRadioReco.framework.channel -import NuRadioReco.framework.trigger -import NuRadioReco.modules.channelSignalReconstructor -signal_reconstructor = NuRadioReco.modules.channelSignalReconstructor.channelSignalReconstructor() - -import uproot -# default uproot readout is awkward arrays, but requires awkward package installed. RNOG data format does not require this. Use numpy instead. -uproot.default_library = "np" - -import numpy as np -from NuRadioReco.utilities import units -import sys -import os -import logging -import time -from scipy import interpolate -import six -from collections import OrderedDict -import astropy.time - - -class readRNOGData: - """ - This is the data reader for RNO-G. Reads RNO-G data from ROOT format using uproot - """ - def __init__(self): - self.logger = logging.getLogger("NuRadioReco.readRNOGdata") - self.__id_current_event = None - self.__t = None - self.__sampling_rate = 3.2 * units.GHz #TODO: 3.2 at the beginning of deployment. Will change to 2.4 GHz after firmware update eventually, but info not yet contained in the .root files. Read out once available. - self._iterator_data = None - self._iterator_header = None - self._data_treename = "waveforms" - self._header_treename = "header" - self.n_events = None - self.input_files = [] - - def begin(self, input_files, input_files_header=None): - - """ - Begin function of the RNO-G reader - - Parameters - ---------- - input_files: list of paths to files containing waveforms - input_files_header: list of paths to files containing header, - if None, headers are expected to be stored in the input_files also - """ - - self.__id_current_event = -1 - self.__t = time.time() - - if isinstance(input_files, six.string_types): - input_files = [input_files] - if isinstance(input_files_header, six.string_types): - input_files_header = [input_files_header] - if input_files_header is None: - input_files_header = input_files - - self.input_files = input_files - self.input_files_header = input_files_header - - self.n_events = 0 - # get the total number of events of all input files - for filename in input_files: - file = uproot.open(filename) - if 'combined' in file: - file = file['combined'] - self.n_events += file[self._data_treename].num_entries - self._set_iterators() - - return self.n_events - - def _set_iterators(self, cut=None): - """ - Set uproot iterators to loop over event trees - - Parameters - ---------- - cut: str - cut string to apply (e.g. for initial event selection based on event_number, ... - e.g. "(event_number==1)" or "(run_number==1)&(event_number<10)" - """ - self.__id_current_event = -1 - - datadict = OrderedDict() - for filename in self.input_files: - if 'combined' in uproot.open(filename): - datadict[filename] = 'combined/' + self._data_treename - else: - datadict[filename] = self._data_treename - - headerdict = OrderedDict() - for filename in self.input_files_header: - if 'combined' in uproot.open(filename): - headerdict[filename] = 'combined/' + self._header_treename - else: - headerdict[filename] = self._header_treename - - # iterator over single events (step 1), for event looping in NuRadioReco dataformat - # may restrict which data to read in the iterator by adding second argument - # read_branches = ['run_number', 'event_number', 'station_number', 'radiant_data[24][2048]'] - self._iterator_data = uproot.iterate(datadict, cut=cut,step_size=1, how=dict, library="np") - self._iterator_header = uproot.iterate(headerdict, cut=cut, step_size=1, how=dict, library="np") - - self.uproot_iterator_data = uproot.iterate(datadict, cut=cut, step_size=1000) - self.uproot_iterator_header = uproot.iterate(headerdict, cut=cut, step_size=1000) - - @register_run() - def run(self, channels=np.arange(24), event_numbers=None, run_numbers=None, cut_string=None): - """ - Run function of the RNOG reader - - Parameters - ---------- - n_channels: int - number of RNOG channels to loop over, default 24 - - event_numbers: None or dict - if dict, use a dict with run number as key and list of event numbers as items - - run_numbers: None or list - list of run numbers to select - Caveat: use only if event_numbers are not set - - cut_string: string - selection string for event pre-selection - Cavieat: use only if event_numbers and run_numbers are not set - """ - - # generate cut string based on passed event_numbers or run_numbers parameters - if not run_numbers is None: - event_cuts = "|".join(["(run_number==%i)" for run_number in run_numbers]) - cut_string = "|".join(event_cuts) - if not event_numbers is None: - event_cuts = [] - for run in event_numbers: - events = event_numbers[run] - for event in events: - event_cuts.append("(run_number==%i)&(event_number==%i)" %(run, event)) - cut_string = "|".join(event_cuts) - self.cut_string = cut_string - - self._set_iterators(cut=self.cut_string) - root_trigger_keys = [ - 'trigger_info.rf_trigger', 'trigger_info.force_trigger', - 'trigger_info.pps_trigger', 'trigger_info.ext_trigger', - 'trigger_info.radiant_trigger', 'trigger_info.lt_trigger', - 'trigger_info.surface_trigger' - ] - self.__t = time.time() - # Note: reading single events is inefficient... - # for event_header, event in zip(self._iterator_header, self._iterator_data): - for event_headers, events in zip(self.uproot_iterator_header, self.uproot_iterator_data): - for event_header, event in zip(event_headers, events): - self.__id_current_event += 1 - #if self.__id_current_event >= self.n_events: - # # all events processed, but iterator should stop before anyways. - # break - if self.__id_current_event % 1000 == 0: - progress = 1. * self.__id_current_event / self.n_events - eta = 0 - if self.__id_current_event > 0: - eta = (time.time() - self.__t) / self.__id_current_event * (self.n_events - self.__id_current_event) / 60. - self.logger.warning("reading in event {}/{} ({:.0f}%) ETA: {:.1f} minutes".format(self.__id_current_event, self.n_events, 100 * progress, eta)) - - run_number = event["run_number"] - evt_number = event["event_number"] - station_id = event_header["station_number"] - self.logger.info("Reading Run: {run_number}, Event {evt_number}, Station {station_id}") - - evt = NuRadioReco.framework.event.Event(run_number, evt_number) - station = NuRadioReco.framework.station.Station(station_id) - #TODO in future: do need to apply calibrations? - - unix_time = event_header["trigger_time"] - event_time = astropy.time.Time(unix_time, format='unix') - - station.set_station_time(event_time) - for trigger_key in root_trigger_keys: - try: - has_triggered = bool(event_header[trigger_key]) - trigger = NuRadioReco.framework.trigger.Trigger(trigger_key.split('.')[-1]) - trigger.set_triggered(has_triggered) - station.set_trigger(trigger) - except ValueError: - pass - - radiant_data = event["radiant_data[24][2048]"] # returns array of n_channels, n_points - # Loop over all requested channels in data - for chan in channels: - channel = NuRadioReco.framework.channel.Channel(chan) - - # Get data from array via graph method - voltage = np.array(radiant_data[chan]) * units.mV - #times = np.arange(len(voltage)) * sampling - - if voltage.shape[0] % 2 != 0: - voltage = voltage[:-1] - - #TODO: need to subtract mean... probably not if running signal reconstructor? - #channel.set_trace(voltage-np.mean(voltage), sampling_rate) - channel.set_trace(voltage, self.__sampling_rate) - station.add_channel(channel) - evt.set_station(station) - # we want to have access to basic signal quantities with implementation from NuRadioReco - #TODO: maybe this should be run in external module? - signal_reconstructor.run(evt, station, None) - yield evt - - def get_events(self): - return self.run() - - def get_n_events(self): - return self.n_events - - def get_filenames(self): - return self.input_files - - def end(self): - pass diff --git a/NuRadioReco/modules/io/rno_g/rnogDataReader.py b/NuRadioReco/modules/io/rno_g/rnogDataReader.py deleted file mode 100644 index cea0be080..000000000 --- a/NuRadioReco/modules/io/rno_g/rnogDataReader.py +++ /dev/null @@ -1,126 +0,0 @@ -import numpy as np -import uproot -import NuRadioReco.framework.event -import NuRadioReco.framework.station -import NuRadioReco.framework.channel -from NuRadioReco.utilities import units -import astropy.time -import glob -import logging -import time -from functools import lru_cache -import six -import NuRadioReco.utilities.metaclasses - -logger = logging.getLogger("RNO-G_IO") -# logger.setLevel(logging.DEBUG) - -# @six.add_metaclass(NuRadioReco.utilities.metaclasses.Singleton) # maybe? -class RNOGDataReader: - - def __init__(self, filenames, *args, **kwargs): - logger.debug("Initializing RNOGDataReader") - self.__filenames = filenames - self.__event_ids = None - self.__sampling_rate = 3.2 * units.GHz #TODO: 3.2 at the beginning of deployment. Will change to 2.4 GHz after firmware update eventually, but info not yet contained in the .root files. Read out once available. - self.__parse_event_ids() - self.__i_events_per_file = np.zeros((len(self.__filenames), 2), dtype=int) - i_event = 0 - for i_file, filename in enumerate(filenames): - file = self.__open_file(filename) - events_in_file = file['waveforms'].num_entries - self.__i_events_per_file[i_file] = [i_event, i_event + events_in_file] - i_event += events_in_file - - self._root_trigger_keys = [ - 'trigger_info.rf_trigger', 'trigger_info.force_trigger', - 'trigger_info.pps_trigger', 'trigger_info.ext_trigger', - 'trigger_info.radiant_trigger', 'trigger_info.lt_trigger', - 'trigger_info.surface_trigger' - ] - - def get_filenames(self): - return self.__filenames - - def get_event_ids(self): - if self.__event_ids is None: - return self.__parse_event_ids() - return self.__event_ids - - def __parse_event_ids(self): - logger.debug('Parsing event ids') - event_ids = np.array([], dtype=int) - run_numbers = np.array([], dtype=int) - for filename in self.__filenames: - file = self.__open_file(filename) - event_ids = np.append(event_ids, file['waveforms']['event_number'].array(library='np').astype(int)) - run_numbers = np.append(run_numbers, file['header']['run_number'].array(library='np').astype(int)) - self.__event_ids = np.array([run_numbers, event_ids]).T - - def __open_file(self, filename): - logger.debug("Opening file {}".format(filename)) - file = uproot.open(filename) - if 'combined' in file: - file = file['combined'] - return file - - def get_n_events(self): - return self.get_event_ids().shape[0] - - # @lru_cache(maxsize=1) # probably not actually relevant outside the data viewer? - def get_event_i(self, i_event): - read_time = time.time() - event = NuRadioReco.framework.event.Event(*self.get_event_ids()[i_event]) - for i_file, filename in enumerate(self.__filenames): - if self.__i_events_per_file[i_file, 0] <= i_event < self.__i_events_per_file[i_file, 1]: - i_event_in_file = i_event - self.__i_events_per_file[i_file, 0] - file = self.__open_file(self.__filenames[i_file]) - station = NuRadioReco.framework.station.Station((file['waveforms']['station_number'].array(library='np', entry_start=i_event_in_file, entry_stop=(i_event_in_file+1))[0])) - # station not set properly in first runs, try from header - if station.get_id() == 0 and 'header' in file: - station = NuRadioReco.framework.station.Station((file['header']['station_number'].array(library='np', entry_start=i_event_in_file, entry_stop=(i_event_in_file+1))[0])) - station.set_is_neutrino() - - if 'header' in file: - unix_time = file['header']['readout_time'].array(library='np', entry_start=i_event_in_file, entry_stop=(i_event_in_file+1))[0] - event_time = astropy.time.Time(unix_time, format='unix') - station.set_station_time(event_time) - ### read in basic trigger data - for trigger_key in self._root_trigger_keys: - try: - has_triggered = bool(file['header'][trigger_key].array(library='np', entry_start=i_event_in_file, entry_stop=(i_event_in_file+1))[0]) - trigger = NuRadioReco.framework.trigger.Trigger(trigger_key.split('.')[-1]) - trigger.set_triggered(has_triggered) - # trigger.set_trigger_time(file['header']['trigger_time']) - station.set_trigger(trigger) - except uproot.exceptions.KeyInFileError: - pass - - waveforms = file['waveforms']['radiant_data[24][2048]'].array(library='np', entry_start=i_event_in_file, entry_stop=(i_event_in_file+1)) - for i_channel in range(waveforms.shape[1]): - channel = NuRadioReco.framework.channel.Channel(i_channel) - channel.set_trace(waveforms[0, i_channel]*units.mV, self.__sampling_rate) - station.add_channel(channel) - event.set_station(station) - logger.debug("Spent {:.0f} ms reading event {}".format((time.time()-read_time) * 1e3, i_event)) - return event - return None - - def get_event(self, event_id): - find_event = np.where((self.get_event_ids()[:,0] == event_id[0]) & (self.get_event_ids()[:,1] == event_id[1]))[0] - if len(find_event) == 0: - return None - elif len(find_event) == 1: - return self.get_event_i(find_event[0]) - else: - raise RuntimeError('There are multiple events with the ID [{}, {}] in the file'.format(event_id[0], event_id[1])) - - def get_events(self): - for ev_i in range(self.get_n_events()): - yield self.get_event_i(ev_i) - - def get_detector(self): - return None - - def get_header(self): - return None diff --git a/NuRadioReco/modules/measured_noise/RNO_G/__init__.py b/NuRadioReco/modules/measured_noise/RNO_G/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/NuRadioReco/modules/measured_noise/RNO_G/noiseImporter.py b/NuRadioReco/modules/measured_noise/RNO_G/noiseImporter.py new file mode 100644 index 000000000..c0c781c8b --- /dev/null +++ b/NuRadioReco/modules/measured_noise/RNO_G/noiseImporter.py @@ -0,0 +1,220 @@ +import numpy as np +import glob +import os +import collections +import time + +from NuRadioReco.modules.io.RNO_G.readRNOGDataMattak import readRNOGData +from NuRadioReco.modules.base.module import register_run +from NuRadioReco.utilities import units + +import logging + + +class noiseImporter: + """ + Imports recorded traces from RNOG stations. + + """ + + + def begin(self, noise_folders, file_pattern="*", + match_station_id=False, station_ids=None, + channel_mapping=None, scramble_noise_file_order=True, + log_level=logging.INFO, random_seed=None, reader_kwargs={}): + """ + + Parameters + ---------- + noise_folders: str or list(str) + Folder(s) containing noise file(s). Search in any subfolder as well. + + file_pattern: str + File pattern used to search for directories, (Default: "*", other examples might be "combined") + + match_station_id: bool + If True, add only noise from stations with the same id. (Default: False) + + station_ids: list(int) + Only add noise from those station ids. If None, use any station. (Default: None) + + channel_mapping: dict or None + option relevant for MC studies of new station designs where we do not + have forced triggers for. The channel_mapping dictionary maps the channel + ids of the MC station to the channel ids of the noise data + Default is None which is 1-to-1 mapping + + scramble_noise_file_order: bool + If True, randomize the order of noise files before reading them. (Default: True) + + log_level: loggging log level + The log level to controll verbosity. (Default: logging.INFO) + + random_seed: int + Seed for the random number generator. (Default: None, no fixed seed). + + reader_kwargs: dict + Optional arguements passed to readRNOGDataMattak + """ + + self.logger = logging.getLogger('NuRadioReco.RNOG.noiseImporter') + self.logger.setLevel(log_level) + self.__random_gen = np.random.Generator(np.random.Philox(random_seed)) + + self._match_station_id = match_station_id + self.__station_ids = station_ids + + self.__channel_mapping = channel_mapping + + self.logger.info(f"\n\tMatch station id: {match_station_id}" + f"\n\tUse noise from only those stations: {station_ids}" + f"\n\tUse the following channel mapping: {channel_mapping}" + f"\n\tRandomize sequence of noise files: {scramble_noise_file_order}") + + if not isinstance(noise_folders, list): + noise_folders = [noise_folders] + + # find all subfolders + noise_files = [] + for noise_folder in noise_folders: + if noise_folder == "": + continue + + noise_files += glob.glob(f"{noise_folder}/**/{file_pattern}root", recursive=True) + self.__noise_folders = np.unique([os.path.dirname(e) for e in noise_files]) + + self.logger.info(f"Found {len(self.__noise_folders)}") + if not len(self.__noise_folders): + self.logger.error("No folders found") + raise FileNotFoundError("No folders found") + + if scramble_noise_file_order: + self.__random_gen.shuffle(self.__noise_folders) + + self._noise_reader = readRNOGData() + + default_reader_kwargs = { + "selectors": [lambda einfo: einfo.triggerType == "FORCE"], + "log_level": log_level, "select_runs": True, "max_trigger_rate": 2 * units.Hz, + "run_types": ["physics"] + } + default_reader_kwargs.update(reader_kwargs) + + self._noise_reader.begin(self.__noise_folders, **default_reader_kwargs) + + # instead of reading all noise events into memory we only get certain information here and read all data in run() + self.logger.info("Get event informations ...") + t0 = time.time() + noise_information = self._noise_reader.get_events_information(keys=["station"]) + self.logger.info(f"... of {len(noise_information)} (selected) events in {time.time() - t0:.2f}s") + + self.__event_index_list = np.array(list(noise_information.keys())) + self.__station_id_list = np.array([ele["station"] for ele in noise_information.values()]) + + self._n_use_event = collections.defaultdict(int) + + + def __get_noise_channel(self, channel_id): + if self.__channel_mapping is None: + return channel_id + else: + return self.__channel_mapping[channel_id] + + + def __draw_noise_event(self, mask): + """ + reader.get_event_by_index can return None when, e.g., the trigger time is inf or the sampling rate 0. + Hence, try again if that happens (should only occur rearly). + + Parameters + ---------- + + mask: np.array(bool) + Mask of which noise events are allowed (e.g. because of matching station ids, ...) + + Returns + ------- + + noise_event: NuRadioReco.framework.event + A event containing noise traces + + i_noise: int + The index of the drawn event + """ + tries = 0 + while tries < 100: + # int(..) necessary because pyroot can not handle np.int64 + i_noise = int(self.__random_gen.choice(self.__event_index_list[mask])) + noise_event = self._noise_reader.get_event_by_index(i_noise) + tries += 1 + if noise_event is not None: + break + + if noise_event is None: + err = "Could not draw a random station which is not None after 100 tries. Stop." + self.logger.error(err) + raise ValueError(err) + + self._n_use_event[i_noise] += 1 + return noise_event, i_noise + + + @register_run() + def run(self, evt, station, det): + + if self._match_station_id: + # select only noise events from simulated station id + station_mask = self.__station_id_list == station.get_id() + if not np.any(station_mask): + raise ValueError(f"No station with id {station.get_id()} in noise data.") + + else: + # select all noise events + station_mask = np.full_like(self.__event_index_list, True) + + noise_event, i_noise = self.__draw_noise_event(station_mask) + + station_id = noise_event.get_station_ids()[0] + noise_station = noise_event.get_station(station_id) + + if self.__station_ids is not None and not station_id in self.__station_ids: + raise ValueError(f"Station id {station_id} not in list of allowed ids: {self.__station_ids}") + + self.logger.debug("Selected noise event {} ({}, run {}, event {})".format( + i_noise, noise_station.get_station_time(), noise_event.get_run_number(), + noise_event.get_id())) + + for channel in station.iter_channels(): + channel_id = channel.get_id() + + trace = channel.get_trace() + noise_channel = noise_station.get_channel(self.__get_noise_channel(channel_id)) + noise_trace = noise_channel.get_trace() + + if len(trace) > 2048: + self.logger.warn("Simulated trace is longer than 2048 bins... trim with :2048") + trace = trace[:2048] + + # sanity checks + if len(trace) != len(noise_trace): + erg_msg = f"Mismatch in trace lenght: Noise has {len(noise_trace)} " + \ + "and simulation has {len(trace)} samples" + self.logger.error(erg_msg) + raise ValueError(erg_msg) + + if channel.get_sampling_rate() != noise_channel.get_sampling_rate(): + erg_msg = "Mismatch in sampling rate: Noise has {} and simulation has {} GHz".format( + noise_channel.get_sampling_rate() / units.GHz, channel.get_sampling_rate() / units.GHz) + self.logger.error(erg_msg) + raise ValueError(erg_msg) + + trace = trace + noise_trace + channel.set_trace(trace, channel.get_sampling_rate()) + + def end(self): + self._noise_reader.end() + n_use = np.array(list(self._n_use_event.values())) + sort = np.flip(np.argsort(n_use)) + self.logger.info("\n\tThe five most used noise events have been used: {}" + .format(", ".join([str(ele) for ele in n_use[sort][:5]]))) + pass diff --git a/NuRadioReco/modules/neutrinoDirectionReconstructor/voltageToEfieldAnalyticConverterForNeutrinos.py b/NuRadioReco/modules/neutrinoDirectionReconstructor/voltageToEfieldAnalyticConverterForNeutrinos.py index e672b51d6..07f7daa19 100644 --- a/NuRadioReco/modules/neutrinoDirectionReconstructor/voltageToEfieldAnalyticConverterForNeutrinos.py +++ b/NuRadioReco/modules/neutrinoDirectionReconstructor/voltageToEfieldAnalyticConverterForNeutrinos.py @@ -187,14 +187,14 @@ def minimizer(params, minimizer=True): order = 5 b, a = signal.butter(order, passband_high[use_channels[iA]], 'bandpass', analog=True) w, h = signal.freqs(b, a, ff[mask]) - f = np.zeros_like(ff, dtype=np.complex) + f = np.zeros_like(ff, dtype=complex) f[mask] = h trace_spectrum *= f order = 10 b, a = signal.butter(order, passband_low[use_channels[iA]], 'bandpass', analog=True) w, h = signal.freqs(b, a, ff[mask]) - f = np.zeros_like(ff, dtype=np.complex) + f = np.zeros_like(ff, dtype=complex) f[mask] = h trace_spectrum *= f @@ -352,8 +352,8 @@ def plotSimulatedVersusReconstructedTraces(nu_zenith, nu_azimuth, shower_energy, travel_distance = np.zeros((n_antennas, maxNumRayTracingSolPerChan)) attenuation = np.zeros((n_antennas, maxNumRayTracingSolPerChan, len(ff))) focusing = np.zeros((n_antennas, maxNumRayTracingSolPerChan, 1)) - reflection_coefficients_theta = np.ones((n_antennas, maxNumRayTracingSolPerChan), dtype=np.complex) - reflection_coefficients_phi = np.ones((n_antennas, maxNumRayTracingSolPerChan), dtype=np.complex) + reflection_coefficients_theta = np.ones((n_antennas, maxNumRayTracingSolPerChan), dtype=complex) + reflection_coefficients_phi = np.ones((n_antennas, maxNumRayTracingSolPerChan), dtype=complex) travel_time_min = float('inf') for iA, position in enumerate(antenna_positions): r = ray.ray_tracing(icemodel, attenuation_model=attenuation_model, n_frequencies_integration=25, n_reflections=n_reflections) diff --git a/NuRadioReco/modules/phasedarray/triggerSimulator.py b/NuRadioReco/modules/phasedarray/triggerSimulator.py index 9fba9b25d..44ac9cc53 100644 --- a/NuRadioReco/modules/phasedarray/triggerSimulator.py +++ b/NuRadioReco/modules/phasedarray/triggerSimulator.py @@ -222,7 +222,7 @@ def power_sum(self, coh_sum, window, step, adc_output='voltage'): num_frames = int(np.floor((len(coh_sum) - window) / step)) if(adc_output == 'voltage'): - coh_sum_squared = (coh_sum * coh_sum).astype(np.float) + coh_sum_squared = (coh_sum * coh_sum).astype(float) elif(adc_output == 'counts'): coh_sum_squared = (coh_sum * coh_sum).astype(int) @@ -230,7 +230,7 @@ def power_sum(self, coh_sum, window, step, adc_output='voltage'): (coh_sum_squared.strides[0] * step, coh_sum_squared.strides[0])) power = np.sum(coh_sum_windowed, axis=1) - return power.astype(np.float) / window, num_frames + return power.astype(float) / window, num_frames def phase_signals(self, traces, beam_rolls): """ diff --git a/NuRadioReco/modules/trigger/highLowThreshold.py b/NuRadioReco/modules/trigger/highLowThreshold.py index c08e9d7c4..cff116833 100644 --- a/NuRadioReco/modules/trigger/highLowThreshold.py +++ b/NuRadioReco/modules/trigger/highLowThreshold.py @@ -32,7 +32,7 @@ def get_high_low_triggers(trace, high_threshold, low_threshold, the bins where the trigger condition is satisfied """ n_bins_coincidence = int(np.round(time_coincidence / dt)) + 1 - c = np.ones(n_bins_coincidence, dtype=np.bool) + c = np.ones(n_bins_coincidence, dtype=bool) logger.debug("length of trace {} bins, coincidence window {} bins".format(len(trace), len(c))) c2 = np.array([1, -1]) @@ -70,7 +70,7 @@ def get_majority_logic(tts, number_of_coincidences=2, time_coincidence=32 * unit if(n_bins_coincidence > n): # reduce coincidence window to maximum trace length n_bins_coincidence = n logger.debug("specified coincidence window longer than tracelenght, reducing coincidence window to trace length") - c = np.ones(n_bins_coincidence, dtype=np.bool) + c = np.ones(n_bins_coincidence, dtype=bool) for i in range(len(tts)): logger.debug("get_majority_logic() length of trace {} bins, coincidence window {} bins".format(len(tts[i]), len(c))) diff --git a/NuRadioReco/modules/voltageToEfieldConverter.py b/NuRadioReco/modules/voltageToEfieldConverter.py index 2b66efa12..df6144a48 100644 --- a/NuRadioReco/modules/voltageToEfieldConverter.py +++ b/NuRadioReco/modules/voltageToEfieldConverter.py @@ -77,7 +77,7 @@ def get_array_of_channels(station, use_channels, det, zenith, azimuth, for iCh, trace in enumerate(traces): V_timedomain[iCh] = trace.get_trace() frequencies = traces[0].get_frequencies() # assumes that all channels have the same sampling rate - V = np.zeros((len(use_channels), len(frequencies)), dtype=np.complex) + V = np.zeros((len(use_channels), len(frequencies)), dtype=complex) for iCh, trace in enumerate(traces): V[iCh] = trace.get_frequency_spectrum() @@ -170,7 +170,7 @@ def run(self, evt, station, det, use_channels=None, use_MC_direction=False, forc E1[mask] = (V[0] * efield_antenna_factor[-1][1] - V[-1] * efield_antenna_factor[0][1])[mask] / denom[mask] E2[mask] = (V[-1] - efield_antenna_factor[-1][0] * E1)[mask] / efield_antenna_factor[-1][1][mask] # solve it in a vectorized way - efield3_f = np.zeros((2, n_frequencies), dtype=np.complex) + efield3_f = np.zeros((2, n_frequencies), dtype=complex) if force_Polarization == 'eTheta': efield3_f[:1, mask] = np.moveaxis(stacked_lstsq(np.moveaxis(efield_antenna_factor[:, 0, mask], 1, 0)[:, :, np.newaxis], np.moveaxis(V[:, mask], 1, 0)), 0, 1) elif force_Polarization == 'ePhi': @@ -178,7 +178,7 @@ def run(self, evt, station, det, use_channels=None, use_MC_direction=False, forc else: efield3_f[:, mask] = np.moveaxis(stacked_lstsq(np.moveaxis(efield_antenna_factor[:, :, mask], 2, 0), np.moveaxis(V[:, mask], 1, 0)), 0, 1) # add eR direction - efield3_f = np.array([np.zeros_like(efield3_f[0], dtype=np.complex), + efield3_f = np.array([np.zeros_like(efield3_f[0], dtype=complex), efield3_f[0], efield3_f[1]]) diff --git a/NuRadioReco/modules/voltageToEfieldConverterPerChannel.py b/NuRadioReco/modules/voltageToEfieldConverterPerChannel.py index d83319e16..bdc79caf8 100644 --- a/NuRadioReco/modules/voltageToEfieldConverterPerChannel.py +++ b/NuRadioReco/modules/voltageToEfieldConverterPerChannel.py @@ -65,7 +65,7 @@ def run(self, evt, station, det, pol=0): trace = channel.get_frequency_spectrum() mask1 = np.abs(efield_antenna_factor[iCh][0]) != 0 mask2 = np.abs(efield_antenna_factor[iCh][1]) != 0 - efield_spectrum = np.zeros((3, len(trace)), dtype=np.complex) + efield_spectrum = np.zeros((3, len(trace)), dtype=complex) efield_spectrum[1][mask1] = (1.0 - pol) ** 2 * trace[mask1] / efield_antenna_factor[iCh][0][mask1] efield_spectrum[2][mask2] = pol ** 2 * trace[mask2] / efield_antenna_factor[iCh][1][mask2] efield.set_frequency_spectrum(efield_spectrum, sampling_rate) diff --git a/NuRadioReco/test/tiny_reconstruction/reference.json b/NuRadioReco/test/tiny_reconstruction/reference.json index 40c6b9def..968d3584b 100644 --- a/NuRadioReco/test/tiny_reconstruction/reference.json +++ b/NuRadioReco/test/tiny_reconstruction/reference.json @@ -1,344 +1 @@ -{ - "0" : { - "station_parameters" : { - "nu_zenith" : null, - "nu_azimuth" : null, - "nu_energy" : null, - "nu_flavor" : null, - "ccnc" : null, - "nu_vertex" : null, - "inelasticity" : null, - "triggered" : null, - "cr_energy" : null, - "cr_zenith" : null, - "cr_azimuth" : null, - "channels_max_amplitude" : 0.08853022471996891, - "zenith" : 0.79, - "azimuth" : 3.96, - "zenith_cr_templatefit" : null, - "zenith_nu_templatefit" : null, - "cr_xcorrelations" : null, - "nu_xcorrelations" : null, - "station_time" : null, - "cr_energy_em" : null, - "nu_inttype" : null, - "chi2_efield_time_direction_fit" : null, - "ndf_efield_time_direction_fit" : null, - "cr_xmax" : null, - "vertex_2D_fit" : null, - "distance_correlations" : null - }, - "sim_station_parameters" : { - "nu_zenith" : null, - "nu_azimuth" : null, - "nu_energy" : null, - "nu_flavor" : null, - "ccnc" : null, - "nu_vertex" : null, - "inelasticity" : null, - "triggered" : null, - "cr_energy" : 1.58489319246E18, - "cr_zenith" : null, - "cr_azimuth" : null, - "channels_max_amplitude" : null, - "zenith" : 0.7853981633974483, - "azimuth" : 3.957853280977215, - "zenith_cr_templatefit" : null, - "zenith_nu_templatefit" : null, - "cr_xcorrelations" : null, - "nu_xcorrelations" : null, - "station_time" : null, - "cr_energy_em" : 1.39187479744259994E18, - "nu_inttype" : null, - "chi2_efield_time_direction_fit" : null, - "ndf_efield_time_direction_fit" : null, - "cr_xmax" : 646.2024663, - "vertex_2D_fit" : null, - "distance_correlations" : null - }, - "channel_parameters" : { - "zenith" : [ null, null, null, null ], - "azimuth" : [ null, null, null, null ], - "maximum_amplitude" : [ 0.0873141653492205, 0.08853022471996891, 0.07885012049626394, 0.07486104950522432 ], - "SNR" : [ { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.004149416647868447, - "peak_amplitude" : 0.004225325303006149, - "Seckel_2_noise" : 5 - }, { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.0061325770929495016, - "peak_amplitude" : 0.006227798651940198, - "Seckel_2_noise" : 5 - }, { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.0038091911965324092, - "peak_amplitude" : 0.0038590946944723325, - "Seckel_2_noise" : 5 - }, { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.003328271897442733, - "peak_amplitude" : 0.0034583171387571744, - "Seckel_2_noise" : 5 - } ], - "maximum_amplitude_envelope" : [ 0.08772020520899446, 0.08892789858406401, 0.07966830072878535, 0.08297084302728261 ], - "P2P_amplitude" : [ 0.1724152251327864, 0.16341917302469117, 0.146681542086196, 0.1458075155635604 ], - "cr_xcorrelations" : [ null, null, null, null ], - "nu_xcorrelations" : [ null, null, null, null ], - "signal_time" : [ -168.2559785714568, -171.65597857145679, -183.2559785714568, -173.85597857145672 ], - "noise_rms" : [ 0.008836023847618458, 0.008955465322914774, 0.008881819339363746, 0.009304834436981291 ], - "signal_regions" : [ null, null, null, null ], - "noise_regions" : [ null, null, null, null ], - "signal_time_offset" : [ null, null, null, null ], - "signal_receiving_zenith" : [ null, null, null, null ], - "signal_ray_type" : [ null, null, null, null ], - "signal_receiving_azimuth" : [ null, null, null, null ] - }, - "electric_field_parameters" : { - "ray_path_type" : [ null, null ], - "polarization_angle" : [ 1.437014716504628, 1.505812376998251 ], - "polarization_angle_expectation" : [ 1.3259385237455839, -1.8156541298442093 ], - "signal_energy_fluence" : [ [ 0.0, 0.2505871320257097, 13.83446329404234 ], [ 0.0, 0.07493732342984019, 17.695470024342303 ] ], - "cr_spectrum_slope" : [ null, -5.071257963761185 ], - "zenith" : [ 0.79, 0.79 ], - "azimuth" : [ 3.96, 3.96 ], - "signal_time" : [ -254.95450631365475, null ], - "nu_vertex_distance" : [ null, null ], - "nu_viewing_angle" : [ null, null ], - "max_amp_antenna" : [ null, null ], - "max_amp_antenna_envelope" : [ null, null ], - "reflection_coefficient_theta" : [ null, null ], - "reflection_coefficient_phi" : [ null, null ], - "cr_spectrum_quadratic_term" : [ null, 3.7278692487503235 ], - "energy_fluence_ratios" : [ null, null ] - } - }, - "1" : { - "station_parameters" : { - "nu_zenith" : null, - "nu_azimuth" : null, - "nu_energy" : null, - "nu_flavor" : null, - "ccnc" : null, - "nu_vertex" : null, - "inelasticity" : null, - "triggered" : null, - "cr_energy" : null, - "cr_zenith" : null, - "cr_azimuth" : null, - "channels_max_amplitude" : 0.33359702926988355, - "zenith" : 0.78, - "azimuth" : 3.95, - "zenith_cr_templatefit" : null, - "zenith_nu_templatefit" : null, - "cr_xcorrelations" : null, - "nu_xcorrelations" : null, - "station_time" : null, - "cr_energy_em" : null, - "nu_inttype" : null, - "chi2_efield_time_direction_fit" : null, - "ndf_efield_time_direction_fit" : null, - "cr_xmax" : null, - "vertex_2D_fit" : null, - "distance_correlations" : null - }, - "sim_station_parameters" : { - "nu_zenith" : null, - "nu_azimuth" : null, - "nu_energy" : null, - "nu_flavor" : null, - "ccnc" : null, - "nu_vertex" : null, - "inelasticity" : null, - "triggered" : null, - "cr_energy" : 1.58489319246E18, - "cr_zenith" : null, - "cr_azimuth" : null, - "channels_max_amplitude" : null, - "zenith" : 0.7853981633974483, - "azimuth" : 3.957853280977215, - "zenith_cr_templatefit" : null, - "zenith_nu_templatefit" : null, - "cr_xcorrelations" : null, - "nu_xcorrelations" : null, - "station_time" : null, - "cr_energy_em" : 1.39187479744259994E18, - "nu_inttype" : null, - "chi2_efield_time_direction_fit" : null, - "ndf_efield_time_direction_fit" : null, - "cr_xmax" : 646.2024663, - "vertex_2D_fit" : null, - "distance_correlations" : null - }, - "channel_parameters" : { - "zenith" : [ null, null, null, null ], - "azimuth" : [ null, null, null, null ], - "maximum_amplitude" : [ 0.2835856575324592, 0.33359702926988355, 0.2791941553632823, 0.3256776186469762 ], - "SNR" : [ { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.00291206901456176, - "peak_amplitude" : 0.002929725052658198, - "Seckel_2_noise" : 5 - }, { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.00340572893988943, - "peak_amplitude" : 0.0034558407562806964, - "Seckel_2_noise" : 5 - }, { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.006845057653623751, - "peak_amplitude" : 0.00701923788399404, - "Seckel_2_noise" : 5 - }, { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.007164164413376706, - "peak_amplitude" : 0.007202262656811785, - "Seckel_2_noise" : 5 - } ], - "maximum_amplitude_envelope" : [ 0.3003272394395557, 0.3399598814912192, 0.29236891303881996, 0.35467083189765575 ], - "P2P_amplitude" : [ 0.5632116726571561, 0.6637288453900938, 0.554719957942302, 0.6479974651404222 ], - "cr_xcorrelations" : [ null, null, null, null ], - "nu_xcorrelations" : [ null, null, null, null ], - "signal_time" : [ 587.7440214285432, 589.9440214285432, 574.9440214285432, 571.3440214285433 ], - "noise_rms" : [ 0.009006474050937649, 0.008732105958313228, 0.007788519158762668, 0.009397032099363559 ], - "signal_regions" : [ null, null, null, null ], - "noise_regions" : [ null, null, null, null ], - "signal_time_offset" : [ null, null, null, null ], - "signal_receiving_zenith" : [ null, null, null, null ], - "signal_ray_type" : [ null, null, null, null ], - "signal_receiving_azimuth" : [ null, null, null, null ] - }, - "electric_field_parameters" : { - "ray_path_type" : [ null, null ], - "polarization_angle" : [ 1.423905497514318, 1.4409019604509443 ], - "polarization_angle_expectation" : [ 1.3232105351636871, -1.818382118426106 ], - "signal_energy_fluence" : [ [ 0.0, 6.607535268712784, 301.8361992200022 ], [ 0.0, 5.437314688159896, 318.63935252848177 ] ], - "cr_spectrum_slope" : [ null, -6.6760948093020485 ], - "zenith" : [ 0.78, 0.78 ], - "azimuth" : [ 3.95, 3.95 ], - "signal_time" : [ 499.91208204952255, null ], - "nu_vertex_distance" : [ null, null ], - "nu_viewing_angle" : [ null, null ], - "max_amp_antenna" : [ null, null ], - "max_amp_antenna_envelope" : [ null, null ], - "reflection_coefficient_theta" : [ null, null ], - "reflection_coefficient_phi" : [ null, null ], - "cr_spectrum_quadratic_term" : [ null, -2.215049003838984 ], - "energy_fluence_ratios" : [ null, null ] - } - }, - "2" : { - "station_parameters" : { - "nu_zenith" : null, - "nu_azimuth" : null, - "nu_energy" : null, - "nu_flavor" : null, - "ccnc" : null, - "nu_vertex" : null, - "inelasticity" : null, - "triggered" : null, - "cr_energy" : null, - "cr_zenith" : null, - "cr_azimuth" : null, - "channels_max_amplitude" : 0.28471716111983614, - "zenith" : 0.78, - "azimuth" : 3.95, - "zenith_cr_templatefit" : null, - "zenith_nu_templatefit" : null, - "cr_xcorrelations" : null, - "nu_xcorrelations" : null, - "station_time" : null, - "cr_energy_em" : null, - "nu_inttype" : null, - "chi2_efield_time_direction_fit" : null, - "ndf_efield_time_direction_fit" : null, - "cr_xmax" : null, - "vertex_2D_fit" : null, - "distance_correlations" : null - }, - "sim_station_parameters" : { - "nu_zenith" : null, - "nu_azimuth" : null, - "nu_energy" : null, - "nu_flavor" : null, - "ccnc" : null, - "nu_vertex" : null, - "inelasticity" : null, - "triggered" : null, - "cr_energy" : 1.58489319246E18, - "cr_zenith" : null, - "cr_azimuth" : null, - "channels_max_amplitude" : null, - "zenith" : 0.7853981633974483, - "azimuth" : 3.957853280977215, - "zenith_cr_templatefit" : null, - "zenith_nu_templatefit" : null, - "cr_xcorrelations" : null, - "nu_xcorrelations" : null, - "station_time" : null, - "cr_energy_em" : 1.39187479744259994E18, - "nu_inttype" : null, - "chi2_efield_time_direction_fit" : null, - "ndf_efield_time_direction_fit" : null, - "cr_xmax" : 646.2024663, - "vertex_2D_fit" : null, - "distance_correlations" : null - }, - "channel_parameters" : { - "zenith" : [ null, null, null, null ], - "azimuth" : [ null, null, null, null ], - "maximum_amplitude" : [ 0.2617580175047184, 0.27712092224880147, 0.2555336277305896, 0.28471716111983614 ], - "SNR" : [ { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.0043472269995871535, - "peak_amplitude" : 0.0045710624556043215, - "Seckel_2_noise" : 5 - }, { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.003710424290754701, - "peak_amplitude" : 0.0037252751132054716, - "Seckel_2_noise" : 5 - }, { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.00777307702513199, - "peak_amplitude" : 0.007907951982053118, - "Seckel_2_noise" : 5 - }, { - "integrated_power" : 0.0, - "peak_2_peak_amplitude" : 0.008133561831805624, - "peak_amplitude" : 0.008336373693178935, - "Seckel_2_noise" : 5 - } ], - "maximum_amplitude_envelope" : [ 0.2653837960611956, 0.2969067177063586, 0.270376775534144, 0.298365691136913 ], - "P2P_amplitude" : [ 0.5188798774102936, 0.550291044285413, 0.5046288232897831, 0.5481116482309814 ], - "cr_xcorrelations" : [ null, null, null, null ], - "nu_xcorrelations" : [ null, null, null, null ], - "signal_time" : [ -134.05597857145676, -132.85597857145672, -146.85597857145672, -144.85597857145672 ], - "noise_rms" : [ 0.009190485861932177, 0.009505810656701827, 0.00886260343524874, 0.0091532389582909 ], - "signal_regions" : [ null, null, null, null ], - "noise_regions" : [ null, null, null, null ], - "signal_time_offset" : [ null, null, null, null ], - "signal_receiving_zenith" : [ null, null, null, null ], - "signal_ray_type" : [ null, null, null, null ], - "signal_receiving_azimuth" : [ null, null, null, null ] - }, - "electric_field_parameters" : { - "ray_path_type" : [ null, null ], - "polarization_angle" : [ 1.4819734258513997, 1.5006361258682528 ], - "polarization_angle_expectation" : [ 1.3232105351636871, -1.818382118426106 ], - "signal_energy_fluence" : [ [ 0.0, 2.5899529944532245, 326.5528775660024 ], [ 0.0, 2.0442071775847332, 413.9200130611196 ] ], - "cr_spectrum_slope" : [ null, -1.7873921532547352 ], - "zenith" : [ 0.78, 0.78 ], - "azimuth" : [ 3.95, 3.95 ], - "signal_time" : [ -193.88791795047757, null ], - "nu_vertex_distance" : [ null, null ], - "nu_viewing_angle" : [ null, null ], - "max_amp_antenna" : [ null, null ], - "max_amp_antenna_envelope" : [ null, null ], - "reflection_coefficient_theta" : [ null, null ], - "reflection_coefficient_phi" : [ null, null ], - "cr_spectrum_quadratic_term" : [ null, 0.4769503271885924 ], - "energy_fluence_ratios" : [ null, null ] - } - } -} \ No newline at end of file +{"0": {"station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": null, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": 0.09219726935932945, "zenith": 0.76585400390625, "azimuth": 3.94125048828125, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": null, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": null, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "sim_station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": 1.58489319246e+18, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": null, "zenith": 0.7853981633974483, "azimuth": 3.957853280977215, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": 1.3918747974426e+18, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": 646.2024663, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "channel_parameters": {"zenith": [null, null, null, null], "azimuth": [null, null, null, null], "maximum_amplitude": [0.0741114148547873, 0.07822416425569524, 0.08019365868532256, 0.09219726935932945], "SNR": [{"integrated_power": 0.0, "peak_2_peak_amplitude": 0.004780092071917627, "peak_amplitude": 0.00485887816779479, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.0029781482057526206, "peak_amplitude": 0.0030899574970230446, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.002963497087732319, "peak_amplitude": 0.0029912274254951183, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.003315202825440897, "peak_amplitude": 0.0035027540130119802, "Seckel_2_noise": 5}], "maximum_amplitude_envelope": [0.0789854151249416, 0.0850534211758587, 0.08035990802734218, 0.09425391321371974], "P2P_amplitude": [0.13819712244096005, 0.15451260952068696, 0.15292399225172826, 0.17650992826729145], "cr_xcorrelations": [null, null, null, null], "nu_xcorrelations": [null, null, null, null], "signal_time": [-165.45597857145674, -169.65597857145679, -182.85597857145672, -181.05597857145676], "noise_rms": [0.009497729224153949, 0.008938569434903116, 0.008083395580028708, 0.009418626437748765], "signal_regions": [null, null, null, null], "noise_regions": [null, null, null, null], "signal_time_offset": [null, null, null, null], "signal_receiving_zenith": [null, null, null, null], "signal_ray_type": [null, null, null, null], "signal_receiving_azimuth": [null, null, null, null]}, "electric_field_parameters": {"ray_path_type": [null, null], "polarization_angle": [1.4654396095601694, 1.4895721505795514], "polarization_angle_expectation": [1.3193420263512765, -1.8222506272385168], "signal_energy_fluence": [[0.0, 0.1499228726724943, 13.406681586865194], [0.0, 0.1138089722425281, 17.17484362157034]], "cr_spectrum_slope": [null, -5.857288373436567], "zenith": [0.76585400390625, 0.76585400390625], "azimuth": [3.94125048828125, 3.94125048828125], "signal_time": [-256.0418188450112, null], "nu_vertex_distance": [null, null], "nu_viewing_angle": [null, null], "max_amp_antenna": [null, null], "max_amp_antenna_envelope": [null, null], "reflection_coefficient_theta": [null, null], "reflection_coefficient_phi": [null, null], "cr_spectrum_quadratic_term": [null, 5.585796124606195], "energy_fluence_ratios": [null, null]}}, "1": {"station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": null, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": 0.34265971445728505, "zenith": 0.79, "azimuth": 3.96, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": null, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": null, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "sim_station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": 1.58489319246e+18, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": null, "zenith": 0.7853981633974483, "azimuth": 3.957853280977215, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": 1.3918747974426e+18, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": 646.2024663, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "channel_parameters": {"zenith": [null, null, null, null], "azimuth": [null, null, null, null], "maximum_amplitude": [0.29720895105426465, 0.34265971445728505, 0.28631950709258547, 0.3355886336552229], "SNR": [{"integrated_power": 0.0, "peak_2_peak_amplitude": 0.006919176246584566, "peak_amplitude": 0.007147260712361315, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.0047795421977155994, "peak_amplitude": 0.00483475221965112, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.007958346991753162, "peak_amplitude": 0.00808818541204719, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.006171263679930979, "peak_amplitude": 0.00634102083568216, "Seckel_2_noise": 5}], "maximum_amplitude_envelope": [0.30136927893948445, 0.3430368419418388, 0.2883085667945778, 0.3430413736731565], "P2P_amplitude": [0.5781561921741288, 0.6591005848543372, 0.5651166372487595, 0.6691483403791307], "cr_xcorrelations": [null, null, null, null], "nu_xcorrelations": [null, null, null, null], "signal_time": [583.9440214285432, 585.5440214285434, 570.5440214285434, 575.1440214285433], "noise_rms": [0.009359269057858015, 0.009003414652873262, 0.008341757854147532, 0.009479378300635246], "signal_regions": [null, null, null, null], "noise_regions": [null, null, null, null], "signal_time_offset": [null, null, null, null], "signal_receiving_zenith": [null, null, null, null], "signal_ray_type": [null, null, null, null], "signal_receiving_azimuth": [null, null, null, null]}, "electric_field_parameters": {"ray_path_type": [null, null], "polarization_angle": [1.4164838522395478, 1.4309598632203113], "polarization_angle_expectation": [1.3259385237455839, -1.8156541298442093], "signal_energy_fluence": [[0.0, 7.388558320124647, 305.36881995591784], [0.0, 6.420446226373949, 324.0685176647595]], "cr_spectrum_slope": [null, -6.704957977720079], "zenith": [0.79, 0.79], "azimuth": [3.96, 3.96], "signal_time": [500.04549368634525, null], "nu_vertex_distance": [null, null], "nu_viewing_angle": [null, null], "max_amp_antenna": [null, null], "max_amp_antenna_envelope": [null, null], "reflection_coefficient_theta": [null, null], "reflection_coefficient_phi": [null, null], "cr_spectrum_quadratic_term": [null, -2.2958549567008015], "energy_fluence_ratios": [null, null]}}, "2": {"station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": null, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": 0.29293556379400426, "zenith": 0.78, "azimuth": 3.95, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": null, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": null, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "sim_station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": 1.58489319246e+18, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": null, "zenith": 0.7853981633974483, "azimuth": 3.957853280977215, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": 1.3918747974426e+18, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": 646.2024663, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "channel_parameters": {"zenith": [null, null, null, null], "azimuth": [null, null, null, null], "maximum_amplitude": [0.2565625777544711, 0.2778654985637024, 0.2672414043783895, 0.29293556379400426], "SNR": [{"integrated_power": 0.0, "peak_2_peak_amplitude": 0.007948759818769018, "peak_amplitude": 0.008209754518559166, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.002601223824751355, "peak_amplitude": 0.0026817225129317146, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.008923287811149969, "peak_amplitude": 0.009007262903406938, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.0068025654072232655, "peak_amplitude": 0.007054720375806741, "Seckel_2_noise": 5}], "maximum_amplitude_envelope": [0.2618281394747253, 0.28224029243565324, 0.26779467480180347, 0.2932426286212969], "P2P_amplitude": [0.5044222411860875, 0.547608219972491, 0.5318703949237792, 0.5851999240317709], "cr_xcorrelations": [null, null, null, null], "nu_xcorrelations": [null, null, null, null], "signal_time": [-102.65597857145667, -133.45597857145674, -148.2559785714567, -147.45597857145674], "noise_rms": [0.008359098917345579, 0.00848077598024109, 0.008989153845587915, 0.008201090181413125], "signal_regions": [null, null, null, null], "noise_regions": [null, null, null, null], "signal_time_offset": [null, null, null, null], "signal_receiving_zenith": [null, null, null, null], "signal_ray_type": [null, null, null, null], "signal_receiving_azimuth": [null, null, null, null]}, "electric_field_parameters": {"ray_path_type": [null, null], "polarization_angle": [1.4786245054415663, 1.5037562254029073], "polarization_angle_expectation": [1.3232105351636871, -1.818382118426106], "signal_energy_fluence": [[0.0, 2.8633662095356187, 335.1320229092165], [0.0, 1.9443851451542769, 431.3306782776532]], "cr_spectrum_slope": [null, -1.715922037957672], "zenith": [0.78, 0.78], "azimuth": [3.95, 3.95], "signal_time": [-193.88791795047757, null], "nu_vertex_distance": [null, null], "nu_viewing_angle": [null, null], "max_amp_antenna": [null, null], "max_amp_antenna_envelope": [null, null], "reflection_coefficient_theta": [null, null], "reflection_coefficient_phi": [null, null], "cr_spectrum_quadratic_term": [null, 0.5209156977104409], "energy_fluence_ratios": [null, null]}}} \ No newline at end of file diff --git a/NuRadioReco/test/trigger_tests/compare_to_reference.py b/NuRadioReco/test/trigger_tests/compare_to_reference.py index c11d604bd..7cfcdbd5e 100644 --- a/NuRadioReco/test/trigger_tests/compare_to_reference.py +++ b/NuRadioReco/test/trigger_tests/compare_to_reference.py @@ -41,7 +41,7 @@ for prop in properties: if(prop == "trigger_time"): try: - np.testing.assert_allclose(np.array(trigger_results[trigger_name][prop], dtype=np.float64), np.array(reference[trigger_name][prop], dtype=np.float64)) + np.testing.assert_allclose(np.array(trigger_results[trigger_name][prop], dtype=float), np.array(reference[trigger_name][prop], dtype=float)) except AssertionError as e: print('Property {} of trigger {} differs from reference'.format(prop, trigger_name)) print(e) diff --git a/NuRadioReco/utilities/bandpass_filter.py b/NuRadioReco/utilities/bandpass_filter.py index d02035e8a..b91700b55 100644 --- a/NuRadioReco/utilities/bandpass_filter.py +++ b/NuRadioReco/utilities/bandpass_filter.py @@ -48,21 +48,21 @@ def get_filter_response(frequencies, passband, filter_type, order, rp=None): f[np.where(frequencies > passband[1])] = 0. return f elif (filter_type == 'butter'): - f = np.zeros_like(frequencies, dtype=np.complex) + f = np.zeros_like(frequencies, dtype=complex) mask = frequencies > 0 b, a = scipy.signal.butter(order, *scipy_args, analog=True) w, h = scipy.signal.freqs(b, a, frequencies[mask]) f[mask] = h return f elif (filter_type == 'butterabs'): - f = np.zeros_like(frequencies, dtype=np.complex) + f = np.zeros_like(frequencies, dtype=complex) mask = frequencies > 0 b, a = scipy.signal.butter(order, *scipy_args, analog=True) w, h = scipy.signal.freqs(b, a, frequencies[mask]) f[mask] = h return np.abs(f) elif(filter_type == 'cheby1'): - f = np.zeros_like(frequencies, dtype=np.complex) + f = np.zeros_like(frequencies, dtype=complex) mask = frequencies > 0 b, a = scipy.signal.cheby1(order, rp, *scipy_args, analog=True) w, h = scipy.signal.freqs(b, a, frequencies[mask]) diff --git a/NuRadioReco/utilities/noise.py b/NuRadioReco/utilities/noise.py index 2dc66b775..0fd45282b 100644 --- a/NuRadioReco/utilities/noise.py +++ b/NuRadioReco/utilities/noise.py @@ -2,9 +2,11 @@ from NuRadioReco.modules import channelGenericNoiseAdder from NuRadioReco.utilities import units, fft from NuRadioReco.modules.trigger.highLowThreshold import get_high_low_triggers -from NuRadioReco.detector import detector +from NuRadioReco.detector import generic_detector as detector from scipy import constants import datetime +import scipy +import scipy.signal import copy import time @@ -31,10 +33,11 @@ def rolled_sum_roll(traces, rolling): # assume first trace always has no rolling sumtr = traces[0].copy() - for i in range(1,len(traces)): + for i in range(1, len(traces)): sumtr += np.roll(traces[i], rolling[i]) return sumtr + def rolling_indices(traces, rolling): """ pre calculates rolling index array for rolled sum via take @@ -53,6 +56,7 @@ def rolling_indices(traces, rolling): rolling_indices.append(np.roll(idx, roll)) return np.array(rolling_indices).astype(int) + def rolled_sum_take(traces, rolling_indices): """ calculates rolled sum via np.take @@ -72,10 +76,11 @@ def rolled_sum_take(traces, rolling_indices): # assume first trace always has no rolling sumtr = traces[0].copy() - for i in range(1,len(traces)): + for i in range(1, len(traces)): sumtr += np.take(traces[i], rolling_indices[i]) return sumtr + def rolled_sum_slicing(traces, rolling): """ calculates rolled sum via slicing @@ -95,7 +100,7 @@ def rolled_sum_slicing(traces, rolling): # assume first trace always has no rolling sumtr = traces[0].copy() - for i in range(1,len(traces)): + for i in range(1, len(traces)): r = rolling[i] if r > 0: sumtr[:-r] += traces[i][r:] @@ -108,8 +113,6 @@ def rolled_sum_slicing(traces, rolling): return sumtr - - class thermalNoiseGenerator(): def __init__(self, n_samples, sampling_rate, Vrms, threshold, time_coincidence, n_majority, time_coincidence_majority, @@ -228,7 +231,8 @@ class thermalNoiseGeneratorPhasedArray(): def __init__(self, detector_filename, station_id, triggered_channels, Vrms, threshold, ref_index, - noise_type="rayleigh"): + noise_type="rayleigh", log_level=logging.WARNING, + pre_trigger_time=100 * units.ns, trace_length=512 * units.ns): """ Efficient algorithms to generate thermal noise fluctuations that fulfill a phased array trigger @@ -252,19 +256,26 @@ def __init__(self, detector_filename, station_id, triggered_channels, the type of the noise, can be * "rayleigh" (default) * "noise" + pre_trigger_time: float + the time in the trace before the trigger happens + trace_length: float + the total trace length """ + logger.setLevel(log_level) self.debug = False self.max_amp = 0 self.upsampling = 2 - self.det = detector.Detector(json_filename=detector_filename) + self.det = detector.GenericDetector(json_filename=detector_filename) self.det.update(datetime.datetime(2018, 10, 1)) self.n_samples = self.det.get_number_of_samples(station_id, triggered_channels[0]) # assuming same settings for all channels self.sampling_rate = self.det.get_sampling_frequency(station_id, triggered_channels[0]) + self.pre_trigger_bins = int(pre_trigger_time * self.sampling_rate) + self.n_samples_trigger = int(trace_length * self.sampling_rate) det_channel = self.det.get_channel(station_id, triggered_channels[0]) - self.adc_n_bits = det_channel["adc_nbits"] - self.adc_noise_n_bits = det_channel["adc_noise_nbits"] + self.adc_n_bits = det_channel["trigger_adc_nbits"] + self.adc_noise_n_bits = det_channel["trigger_adc_noise_nbits"] self.n_channels = len(triggered_channels) self.triggered_channels = triggered_channels @@ -293,7 +304,6 @@ def __init__(self, detector_filename, station_id, triggered_channels, roll = np.array(np.round(np.array(delays) * self.sampling_rate * self.upsampling)).astype(int) self.beam_time_delays[iBeam] = roll - print(self.beam_time_delays) self.Vrms = Vrms self.threshold = threshold self.noise_type = noise_type @@ -307,10 +317,10 @@ def __init__(self, detector_filename, station_id, triggered_channels, self.filt = channelBandPassFilter.get_filter(self.ff, station_id, channel_id, self.det, passband=[96 * units.MHz, 100 * units.GHz], filter_type='cheby1', order=4, rp=0.1) self.filt *= channelBandPassFilter.get_filter(self.ff, station_id, channel_id, self.det, - passband=[0 * units.MHz, 220 * units.MHz], filter_type='cheby1', order=7, rp=0.1) + passband=[1 * units.MHz, 220 * units.MHz], filter_type='cheby1', order=7, rp=0.1) self.norm = np.trapz(np.abs(self.filt) ** 2, self.ff) self.amplitude = (self.max_freq - self.min_freq) ** 0.5 / self.norm ** 0.5 * self.Vrms - print(f"Vrms = {self.Vrms:.2f}, noise amplitude = {self.amplitude:.2f}, bandwidth = {self.norm/units.MHz:.0f}MHz") + print(f"Vrms = {self.Vrms:.3g}V, noise amplitude = {self.amplitude:.3g}V, bandwidth = {self.norm/units.MHz:.0f}MHz") print(f"frequency range {self.min_freq/units.MHz}MHz - {self.max_freq/units.MHz}MHz") self.adc_ref_voltage = self.Vrms * (2 ** (self.adc_n_bits - 1) - 1) / (2 ** (self.adc_noise_n_bits - 1) - 1) @@ -325,15 +335,14 @@ def __init__(self, detector_filename, station_id, triggered_channels, self.sampling_rate * self.upsampling, self.amplitude, self.noise_type) - def __generation(self): """ separated trace generation part for PA noise trigger """ for iCh in range(self.n_channels): - #spec = self.noise.bandlimited_noise(self.min_freq, self.max_freq, self.n_samples * self.upsampling, + # spec = self.noise.bandlimited_noise(self.min_freq, self.max_freq, self.n_samples * self.upsampling, # self.sampling_rate * self.upsampling, # self.amplitude, self.noise_type, time_domain=False) - + # function that does not re-calculate parameters in each simulated trace spec = self.noise.bandlimited_noise_from_precalculated_parameters(self.noise_type, time_domain=False) spec *= self.filt @@ -341,7 +350,6 @@ def __generation(self): self._traces[iCh] = perfect_floor_comparator(trace, self.adc_n_bits, self.adc_ref_voltage) - def __phasing(self): """ separated phasing part for PA noise trigger """ @@ -362,33 +370,37 @@ def __triggering(self): self.max_amp = 0 # take square over entire array - coh_sum_squared = self._phased_traces**2 + coh_sum_squared = self._phased_traces ** 2 # bin the data into windows of length self.step and normalise to step length - reduced_array = np.add.reduceat(coh_sum_squared.T,np.arange(0,np.shape(coh_sum_squared)[1],self.step)).T / self.step + reduced_array = np.add.reduceat(coh_sum_squared.T, np.arange(0, np.shape(coh_sum_squared)[1], self.step)).T / self.step sliding_windows = [] # self.window can extend over multiple steps, # assuming self.window being an integer multiple of self.step the reduction sums over subsequent steps - steps_per_window = self.window//self.step + steps_per_window = self.window // self.step # better extend the array in order to also trigger on sum of last/first (matching a previous implementation) - extended_reduced_array = np.column_stack([reduced_array, reduced_array[:,0:steps_per_window]]) + extended_reduced_array = np.column_stack([reduced_array, reduced_array[:, 0:steps_per_window]]) for offset in range(steps_per_window): # sum over steps_per_window adjacent steps window_sum = np.add.reduceat(extended_reduced_array.T, np.arange(offset, np.shape(extended_reduced_array)[1], steps_per_window)).T / steps_per_window sliding_windows.append(window_sum) - #self.max_amp = max(np.array(sliding_windows).max(), self.max_amp) + # self.max_amp = max(np.array(sliding_windows).max(), self.max_amp) self.max_amp = np.array(sliding_windows).max() # check if trigger condition is fulfilled anywhere if self.max_amp > self.threshold: - # check in which beam the trigger condition was fulfilled - sliding_windows = np.concatenate(sliding_windows, axis=1) - triggered_beams = np.amax(sliding_windows, axis=1) > self.threshold - for iBeam, is_triggered in enumerate(triggered_beams): - # print out each beam that has triggered - if is_triggered: - logger.info(f"triggered at beam {iBeam}") + sliding_windows = np.array(sliding_windows) + tmp = np.argwhere(sliding_windows > self.threshold) + triggered_step = tmp[0][0] + triggered_bin = tmp[0][2] * (self.window + triggered_step * steps_per_window) if(self.debug): + # check in which beam the trigger condition was fulfilled + sliding_windows = np.concatenate(sliding_windows, axis=1) + triggered_beams = np.amax(sliding_windows, axis=1) > self.threshold + for iBeam, is_triggered in enumerate(triggered_beams): + # print out each beam that has triggered + if is_triggered: + logger.info(f"triggered at beam {iBeam}") import matplotlib.pyplot as plt fig, ax = plt.subplots(5, 1, sharex=True) for iCh in range(self.n_channels): @@ -397,8 +409,8 @@ def __triggering(self): ax[4].plot(self._phased_traces[iBeam]) fig.tight_layout() plt.show() - return True - return False + return True, triggered_bin + return False, None def __triggering_strided(self): """ separated trigger part for PA noise trigger using np.lib.stride_tricks.as_strided """ @@ -411,7 +423,7 @@ def __triggering_strided(self): coh_sum_windowed = np.lib.stride_tricks.as_strided(coh_sum_squared, (num_frames, self.window), (coh_sum_squared.strides[0] * self.step, coh_sum_squared.strides[0])) squared_mean = np.sum(coh_sum_windowed, axis=1) / self.window - #self.max_amp = max(squared_mean.max(), self.max_amp) + # self.max_amp = max(squared_mean.max(), self.max_amp) self.max_amp = max(squared_mean.max(), self.max_amp) if True in (squared_mean > self.threshold): logger.info(f"triggered at beam {iBeam}") @@ -427,7 +439,6 @@ def __triggering_strided(self): return True return False - def generate_noise(self, phasing_mode="slice", trigger_mode="binned_sum", debug=False): """ generates noise traces for all channels that will cause a high/low majority logic trigger @@ -460,7 +471,7 @@ def generate_noise(self, phasing_mode="slice", trigger_mode="binned_sum", debug= while True: counter += 1 if(counter % 1000 == 0): - logger.info(f"{counter:d}, {self.max_amp:.2f}, threshold = {self.threshold:.2f}") + logger.info(f"{counter:d}, {self.max_amp:.2g}, threshold = {self.threshold:.2g}") # some printout for profiling logger.info(f"Time consumption: GENERATION: {dt_generation:.4f}, PHASING: {dt_phasing:.4f}, TRIGGER: {dt_triggering:.4f}") tstart = time.process_time() @@ -480,12 +491,12 @@ def generate_noise(self, phasing_mode="slice", trigger_mode="binned_sum", debug= logger.error(f"Requested phasing_mode {phasing_mode}. Only 'slice' and 'roll' are allowed") raise NotImplementedError(f"Requested phasing_mode {phasing_mode}. Only 'slice' and 'roll' are allowed") - # time profiling phasing + # time profiling phasing dt_phasing += time.process_time() - tstart tstart = time.process_time() if trigger_mode == "binned_sum": - is_triggered = self.__triggering() + is_triggered, triggered_bin = self.__triggering() elif trigger_mode == "stride": # more time consuming attempt to do triggering compared to taking binned sums is_triggered = self.__triggering_strided() @@ -497,7 +508,14 @@ def generate_noise(self, phasing_mode="slice", trigger_mode="binned_sum", debug= dt_triggering += time.process_time() - tstart if is_triggered: - return self._traces, self._phased_traces + triggered_bin = triggered_bin // 2 # the trace is cut in the downsampled version. Therefore, triggered bin is factor of two smaller. + i_low = triggered_bin - self.pre_trigger_bins + i_high = i_low + self.n_samples_trigger + if (i_low >= 0) and (i_high < self.n_samples): + # traces need to be downsampled + # resample and use axis -1 since trace might be either shape (N) for analytic trace or shape (3,N) for E-field + self._traces = scipy.signal.resample(self._traces, np.shape(self._traces)[-1] // self.upsampling, axis=-1) + return self._traces[:, i_low:i_high], self._phased_traces def generate_noise2(self, debug=False): """ @@ -512,7 +530,7 @@ def generate_noise2(self, debug=False): while True: counter += 1 if(counter % 1000 == 0): - print(f"{counter:d}, {max_amp:.2f}, threshold = {self.threshold:.2f}") + print(f"{counter:d}, {max_amp:.3g}, threshold = {self.threshold:.3g}") for iCh in range(self.n_channels): spec = self.noise.bandlimited_noise(self.min_freq, self.max_freq, self.n_samples * self.upsampling, self.sampling_rate * self.upsampling, diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index ca0fe352d..57bc82fb4 100644 --- a/NuRadioReco/utilities/trace_utilities.py +++ b/NuRadioReco/utilities/trace_utilities.py @@ -32,7 +32,7 @@ def get_efield_antenna_factor(station, frequencies, channels, detector, zenith, antenna_pattern_provider: AntennaPatternProvider """ n_ice = ice.get_refractive_index(-0.01, detector.get_site(station.get_id())) - efield_antenna_factor = np.zeros((len(channels), 2, len(frequencies)), dtype=np.complex) # from antenna model in e_theta, e_phi + efield_antenna_factor = np.zeros((len(channels), 2, len(frequencies)), dtype=complex) # from antenna model in e_theta, e_phi for iCh, channel_id in enumerate(channels): zenith_antenna = zenith t_theta = 1. @@ -88,12 +88,12 @@ def get_channel_voltage_from_efield(station, electric_field, channels, detector, spectrum = electric_field.get_frequency_spectrum() efield_antenna_factor = get_efield_antenna_factor(station, frequencies, channels, detector, zenith, azimuth, antenna_pattern_provider) if return_spectrum: - voltage_spectrum = np.zeros((len(channels), len(frequencies)), dtype=np.complex) + voltage_spectrum = np.zeros((len(channels), len(frequencies)), dtype=complex) for i_ch, ch in enumerate(channels): voltage_spectrum[i_ch] = np.sum(efield_antenna_factor[i_ch] * np.array([spectrum[1], spectrum[2]]), axis=0) return voltage_spectrum else: - voltage_trace = np.zeros((len(channels), 2 * (len(frequencies) - 1)), dtype=np.complex) + voltage_trace = np.zeros((len(channels), 2 * (len(frequencies) - 1)), dtype=complex) for i_ch, ch in enumerate(channels): voltage_trace[i_ch] = fft.freq2time(np.sum(efield_antenna_factor[i_ch] * np.array([spectrum[1], spectrum[2]]), axis=0), electric_field.get_sampling_rate()) return np.real(voltage_trace) @@ -218,7 +218,7 @@ def apply_butterworth(spectrum, frequencies, passband, order=8): The filtered spectrum """ - f = np.zeros_like(frequencies, dtype=np.complex) + f = np.zeros_like(frequencies, dtype=complex) mask = frequencies > 0 b, a = scipy.signal.butter(order, passband, 'bandpass', analog=True) w, h = scipy.signal.freqs(b, a, frequencies[mask]) diff --git a/README.md b/README.md index 799965a22..98bd87634 100644 --- a/README.md +++ b/README.md @@ -32,14 +32,16 @@ NuRadioMC is continuously improved and new features are being added. The followi If you would like to contribute, please contact @cg-laser or @anelles for permissions to work on NuRadioMC. We work with pull requests only that can be merged after review. Also please visit https://nu-radio.github.io/NuRadioMC/Introduction/pages/contributing.html for details on our workflow and coding conventions. + +## Publications builing up on NuRadioMC/Reco NuRadioMC is used in an increasing number of studies. To get an overview for what NuRadioMC can be used for, please have a look at the following publications or see [here](https://inspirehep.net/literature?sort=mostrecent&size=25&page=1&q=refersto%3Arecid%3A1738571%20or%20refersto%3Arecid%3A1725583): -* V. B. Valera, M. Bustamante and C. Glaser, “Near-future discovery of the diffuse flux of ultra-high-energy cosmic neutrinos”, [arXiv:2210.03756](https://arxiv.org/abs/2210.03756) -* Alfonso Garcia Soto, Diksha Garg, Mary Hall Reno, Carlos A. Argüelles, "Probing Quantum Gravity with Elastic Interactions of Ultra-High-Energy Neutrinos", [arXiv:2209.06282](https://arxiv.org/abs/2209.06282) -* Damiano F. G. Fiorillo, Mauricio Bustamante, Victor B. Valera, "Near-future discovery of point sources of ultra-high-energy neutrinos", [arXiv:2205.15985](https://arxiv.org/abs/2205.15985) +* V. B. Valera, M. Bustamante and C. Glaser, “Near-future discovery of the diffuse flux of ultra-high-energy cosmic neutrinos”, Phys. Rev. D 107, 043019 [arXiv:2210.03756](https://arxiv.org/abs/2210.03756) +* Alfonso Garcia Soto, Diksha Garg, Mary Hall Reno, Carlos A. Argüelles, "Probing Quantum Gravity with Elastic Interactions of Ultra-High-Energy Neutrinos", Phys. Rev. D 107, 033009 (2023) [arXiv:2209.06282](https://arxiv.org/abs/2209.06282) +* Damiano F. G. Fiorillo, Mauricio Bustamante, Victor B. Valera, "Near-future discovery of point sources of ultra-high-energy neutrinos", JCAP 03 (2023) 026 [arXiv:2205.15985](https://arxiv.org/abs/2205.15985) * C. Glaser, S. McAleer, S. Stjärnholm, P. Baldi, S. W. Barwick, “Deep learning reconstruction of the neutrino direction and energy from in-ice radio detector data”, Astroparticle Physics 145, (2023) 102781, [doi:10.1016/j.astropartphys.2022.102781](https://doi.org/10.1016/j.astropartphys.2022.102781), [arXiv:2205.15872](https://arxiv.org/abs/2205.15872) -* J. Beise and C. Glaser, “In-situ calibration system for the measurement of the snow accumulation and the index-of-refraction profile for radio neutrino detectors”, [arXiv:2205.00726](https://arxiv.org/abs/2205.00726) +* J. Beise and C. Glaser, “In-situ calibration system for the measurement of the snow accumulation and the index-of-refraction profile for radio neutrino detectors”, Journal of Instrumentation 18 P01036 (2023), [arXiv:2205.00726](https://arxiv.org/abs/2205.00726) * V. B. Valera, M. Bustamante and C. Glaser, “The ultra-high-energy neutrino-nucleon cross section: measurement forecasts for an era of cosmic EeV-neutrino discovery”, Journal of High Energy Physics 06 (2022) 105, [doi:10.1007/JHEP06(2022)105](https://doi.org/10.1007/JHEP06(2022%29105), [arXiv:2204.04237](https://arxiv.org/abs/2204.04237) * ARIANNA collaboration (A. Anker et al.), “Measuring the Polarization Reconstruction Resolution of the ARIANNA Neutrino Detector with Cosmic Rays”, Journal of Cosmology and Astroparticle Physics 04(2022)022, [doi:10.1088/1475-7516/2022/04/022](https://doi.org/10.1088/1475-7516/2022/04/022), [arXiv:2112.01501](https://arxiv.org/abs/2112.01501) * ARIANNA collaboration (A. Anker et al.), “Improving sensitivity of the ARIANNA detector by rejecting thermal noise with deep learning”, Journal of Instrumentation 17 P03007 (2022), [doi:10.1088/1748-0221/17/03/P03007](https://doi.org/10.1088/1748-0221/17/03/P03007), [arXiv:2112.01031](https://arxiv.org/abs/2112.01031) diff --git a/changelog.txt b/changelog.txt index 1004f2da1..b30a75ed7 100644 --- a/changelog.txt +++ b/changelog.txt @@ -4,17 +4,25 @@ please update the categories "new features" and "bugfixes" before a pull request version 2.2.0-dev new features: +- added getting to query ARZ charge-excess profiles - upgrade to proposal v 7.5.0 with new NuRadioProposal interface and improved LPM treatment - add script to generate / download pre-calculated proposal tables for known detector configurations - adding default RadioPropa ice model object to medium with the feature to set a personlised object as alternative - add positions array functionality to simple ice model in average and gradient functions - analytic ray tracing solutions are now sorted consistently from lowest to highest ray - added ability to generate high-low-triggered noise on a narrow band but return full-band waveforms +- phased array noise generation utility class: calculate trigger time and cut trace accordingly +- use Philox noise generator in noise adder module (this changes the default random sequence) - added 'block offset' removal/simulation module for RNO-G + bugfixes: - fixed/improved C++ raytracer not finding solutions for some near-horizontal or near-shadowzone vertices - fixed wrong number in Feldman-Cousins upper limit +version 2.1.8 +bugfixes: +- replace deprecated np.float with float + version 2.1.7 new features: - add attenuation model from the 2021 measurements taken at Summit Station diff --git a/pyproject.toml b/pyproject.toml index d936b0ee8..e5e0b4292 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ tinydb = ">=4.1.1" tinydb-serialization = ">=2.1" aenum = "*" astropy = "*" -radiotools = ">=0.2" +radiotools = ">=0.2.1" cython = "*" dash = ">=2.0" future = "*" @@ -50,9 +50,13 @@ proposal = "7.5.1" pygdsm = {git = "https://github.com/telegraphic/pygdsm"} nifty5 = {git = "https://gitlab.mpcdf.mpg.de/ift/nifty.git", branch="NIFTy_5"} pypocketfft = {git = "https://gitlab.mpcdf.mpg.de/mtr/pypocketfft"} +pandas = "*" +mattak = {git = "https://github.com/RNO-G/mattak"} +runtable = {git = "ssh://git@github.com/RNO-G/rnog-runtable.git"} [tool.poetry.extras] documentation = ["Sphinx", "sphinx-rtd-theme", "numpydoc"] proposal = ["proposal"] galacticnoise = ['pygdsm'] ift_reco = ['nifty5', 'pypocketfft'] +RNO_G_DATA = ["mattak", "runtable", "pandas"]