diff --git a/NuRadioReco/modules/io/coreas/coreas.py b/NuRadioReco/modules/io/coreas/coreas.py index 6752490a2..c82490f2f 100644 --- a/NuRadioReco/modules/io/coreas/coreas.py +++ b/NuRadioReco/modules/io/coreas/coreas.py @@ -553,6 +553,12 @@ def initialize_star_shape(self): logger.info(f'Initialize star shape pattern for interpolation. The shower arrives at zenith={self.zenith/units.deg:.0f}deg, azimuth={self.azimuth/units.deg:.0f}deg with radius {self.star_radius:.0f}m in the shower plane and {self.geo_star_radius:.0f}m on ground. ') self.star_shape_initilized = True + def get_sampling_rate(self): + """ + returns the sampling rate of the electric field + """ + return self.sampling_rate + def get_empty_efield(self): """ returns the an array of zeros in the shape of the electric field on the sky @@ -728,11 +734,13 @@ def get_interp_efield_value(self, position_on_ground, core): efield_interp : float interpolated efield value or efield of clostest observer """ + logger.info(f"get interpolated efield for antenna position {position_on_ground} on ground and core position {core}") antenna_position = copy.copy(position_on_ground) #core and antenna need to be in the same z plane antenna_position[2] = core[2] # transform antenna position into shower plane with respect to core position, core position is set to 0,0 in shower plane antenna_pos_vBvvB = self.cs.transform_to_vxB_vxvxB(antenna_position, core=core) + logger.info(f"antenna position in shower plane {antenna_pos_vBvvB}") # calculate distance between core position at(0,0) and antenna positions in shower plane dcore_vBvvB = np.linalg.norm(antenna_pos_vBvvB[:-1]) diff --git a/NuRadioReco/modules/io/coreas/readCoREASDetector.py b/NuRadioReco/modules/io/coreas/readCoREASDetector.py index 19c122717..2d9ab3e26 100644 --- a/NuRadioReco/modules/io/coreas/readCoREASDetector.py +++ b/NuRadioReco/modules/io/coreas/readCoREASDetector.py @@ -14,6 +14,7 @@ import NuRadioReco.framework.radio_shower from NuRadioReco.framework.parameters import showerParameters as shp from NuRadioReco.framework.parameters import stationParameters as stnp +from NuRadioReco.framework.parameters import electricFieldParameters as efp import NuRadioReco.modules.io.coreas from NuRadioReco.modules.io.coreas import coreas from NuRadioReco.utilities import units @@ -32,19 +33,6 @@ def get_random_core_positions(xmin, xmax, ymin, ymax, n_cores, seed=None): np.zeros(n_cores)]).T return cores -def get_efield_times(efield, sampling_rate): - """ - calculate the time axis of the electric field from the sampling rate - - Parameters - ---------- - efield: array (n_samples, n_polarizations) - """ - if efield is None: - return None - efield_times = np.arange(0, len(efield[:,0])) / sampling_rate - return efield_times - def apply_hanning(efields): """ Apply a half hann window to the electric field in the time domain @@ -62,9 +50,9 @@ def apply_hanning(efields): return None else: smoothed_trace = np.zeros_like(efields) - half_hann_window = half_hann_window(efields.shape[0], half_percent=0.1) + filter = half_hann_window(efields.shape[0], half_percent=0.1) for pol in range(efields.shape[1]): - smoothed_trace[:,pol] = efields[:,pol] * half_hann_window + smoothed_trace[:,pol] = efields[:,pol] * filter return smoothed_trace def select_channels_per_station(det, station_id, requested_channel_ids): @@ -102,11 +90,9 @@ def __init__(self): self.__t = 0 self.__t_event_structure = 0 self.__t_per_event = 0 - self.__input_file = None self.__corsika = None self.__interp_lowfreq = None self.__interp_highfreq = None - self.__sampling_rate = None self.logger = logging.getLogger('NuRadioReco.readCoREASDetector') def begin(self, input_file, interp_lowfreq=30*units.MHz, interp_highfreq=1000*units.MHz, log_level=logging.INFO, debug=False): @@ -133,7 +119,7 @@ def begin(self, input_file, interp_lowfreq=30*units.MHz, interp_highfreq=1000*un self.__interp_highfreq = interp_highfreq self.debug = debug - self.coreasInterpolator = NuRadioReco.modules.io.coreas.coreasInterpolator(self.__corsika) + self.coreasInterpolator = coreas.coreasInterpolator(self.__corsika) self.coreasInterpolator.initialize_efield_interpolator(self.__interp_lowfreq, self.__interp_highfreq) @register_run() @@ -158,7 +144,7 @@ def run(self, detector, core_position_list=[], selected_station_ids=[], selected for iCore, core in enumerate(core_position_list): t = time.time() - evt = NuRadioReco.framework.event.Event(evt.get_run_number(), iCore) # create empty event + evt = NuRadioReco.framework.event.Event(self.__corsika.get_run_number(), iCore) # create empty event sim_shower = self.__corsika.get_first_sim_shower() sim_shower.set_parameter(shp.core, core) evt.add_sim_shower(sim_shower) @@ -179,18 +165,23 @@ def run(self, detector, core_position_list=[], selected_station_ids=[], selected selected_channel_ids = channel_ids_in_station channel_ids_dict = select_channels_per_station(detector, station_id, selected_channel_ids) for ch_g_ids in channel_ids_dict.keys(): - antenna_position_rel = detector.get_relative_position(station_id, ch_g_ids) + channel_ids_for_group_id = channel_ids_dict[ch_g_ids] + antenna_position_rel = detector.get_relative_position(station_id, channel_ids_for_group_id[0]) antenna_position = det_station_position + antenna_position_rel res_efield = self.coreasInterpolator.get_interp_efield_value(antenna_position, core) smooth_res_efield = apply_hanning(res_efield) if smooth_res_efield is None: smooth_res_efield = self.coreasInterpolator.get_empty_efield() - efield_times = get_efield_times(smooth_res_efield, self.__sampling_rate) - channel_ids_for_group_id = channel_ids_dict[ch_g_ids] - coreas.add_electric_field_to_sim_station(sim_station, channel_ids_for_group_id, smooth_res_efield.T, efield_times, self.zenith, self.azimuth, self.magnetic_field_vector, self.__sampling_rate ) + electric_field = NuRadioReco.framework.electric_field.ElectricField(channel_ids_for_group_id) + electric_field.set_trace(smooth_res_efield.T, self.coreasInterpolator.get_sampling_rate()) + electric_field.set_trace_start_time(0) # TODO: set this correctly. + electric_field.set_parameter(efp.ray_path_type, 'direct') + electric_field.set_parameter(efp.zenith, sim_shower[shp.zenith]) + electric_field.set_parameter(efp.azimuth, sim_shower[shp.azimuth]) + sim_station.add_electric_field(electric_field) station.set_sim_station(sim_station) - distance_to_core = np.linalg.norm(det_station_position[:-1] - core[:-1]) - station.set_parameter(stnp.distance_to_core, distance_to_core) + # distance_to_core = np.linalg.norm(det_station_position[:-1] - core[:-1]) + # station.set_parameter(stnp.distance_to_core, distance_to_core) evt.set_station(station) self.__t += time.time() - t