diff --git a/NuRadioMC/simulation/simulation.py b/NuRadioMC/simulation/simulation.py index a48a18850..4bddedfe6 100644 --- a/NuRadioMC/simulation/simulation.py +++ b/NuRadioMC/simulation/simulation.py @@ -265,8 +265,9 @@ def __init__(self, inputfilename, ################################ # perfom a dummy detector simulation to determine how the signals are filtered - self._bandwidth_per_channel = {} - self._amplification_per_channel = {} + self._integrated_channel_response = {} + self._integrated_channel_response_normalization = {} + self._max_amplification_per_channel = {} self.__noise_adder_normalization = {} # first create dummy event and station with channels @@ -301,9 +302,13 @@ def __init__(self, inputfilename, self._station.set_station_time(self._evt_time) self._evt.set_station(self._station) + # run detector simulation self._detector_simulation_filter_amp(self._evt, self._station, self._det) - self._bandwidth_per_channel[self._station_id] = {} - self._amplification_per_channel[self._station_id] = {} + + self._integrated_channel_response[self._station_id] = {} + self._integrated_channel_response_normalization[self._station_id] = {} + self._max_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=complex) @@ -311,56 +316,102 @@ def __init__(self, inputfilename, if hasattr(instance, "get_filter"): filt *= instance.get_filter(ff, self._station_id, channel_id, self._det, **kwargs) - self._amplification_per_channel[self._station_id][channel_id] = np.abs(filt).max() - bandwidth = np.trapz(np.abs(filt) ** 2, ff) - self._bandwidth_per_channel[self._station_id][channel_id] = bandwidth - logger.status(f"bandwidth of station {self._station_id} channel {channel_id} is {bandwidth/units.MHz:.1f}MHz") + self._max_amplification_per_channel[self._station_id][channel_id] = np.abs(filt).max() + + mean_integrated_response = np.mean( + np.abs(filt)[np.abs(filt) > np.abs(filt).max() / 100] ** 2) # a factor of 100 corresponds to -40 dB in amplitude + self._integrated_channel_response_normalization[self._station_id][channel_id] = mean_integrated_response + + integrated_channel_response = np.trapz(np.abs(filt) ** 2, ff) + self._integrated_channel_response[self._station_id][channel_id] = integrated_channel_response + + logger.debug(f"Station.channel {self._station_id}.{channel_id} estimated bandwidth is " + f"{integrated_channel_response / mean_integrated_response / units.MHz:.1f} MHz") ################################ - self._bandwidth = next(iter(next(iter(self._bandwidth_per_channel.values())).values())) - amplification = next(iter(next(iter(self._amplification_per_channel.values())).values())) + self._bandwidth = next(iter(next(iter(self._integrated_channel_response.values())).values())) # get value of first station/channel key pair + norm = next(iter(next(iter(self._integrated_channel_response_normalization.values())).values())) + amplification = next(iter(next(iter(self._max_amplification_per_channel.values())).values())) + noise_temp = self._cfg['trigger']['noise_temperature'] Vrms = self._cfg['trigger']['Vrms'] + if noise_temp is not None and Vrms is not None: - raise AttributeError(f"Specifying noise temperature (set to {noise_temp}) and Vrms (set to {Vrms} is not allowed.") + raise AttributeError(f"Specifying noise temperature (set to {noise_temp}) and Vrms (set to {Vrms}) is not allowed).") + + self._Vrms_per_channel = collections.defaultdict(dict) + self._Vrms_efield_per_channel = collections.defaultdict(dict) + if noise_temp is not None: + if noise_temp == "detector": + logger.status("Use noise temperature from detector description to determine noise Vrms in each channel.") self._noise_temp = None # the noise temperature is defined in the detector description else: self._noise_temp = float(noise_temp) - self._Vrms_per_channel = {} - self._noiseless_channels = {} - for station_id in self._bandwidth_per_channel: - self._Vrms_per_channel[station_id] = {} - self._noiseless_channels[station_id] = [] - for channel_id in self._bandwidth_per_channel[station_id]: + logger.status(f"Use a noise temperature of {noise_temp / units.kelvin:.1f} K for each channel to determine noise Vrms.") + + self._noiseless_channels = collections.defaultdict(list) + for station_id in self._integrated_channel_response: + for channel_id in self._integrated_channel_response[station_id]: if self._noise_temp is None: noise_temp_channel = self._det.get_noise_temperature(station_id, channel_id) else: noise_temp_channel = self._noise_temp + if self._det.is_channel_noiseless(station_id, channel_id): self._noiseless_channels[station_id].append(channel_id) - self._Vrms_per_channel[station_id][channel_id] = (noise_temp_channel * 50 * constants.k * - self._bandwidth_per_channel[station_id][channel_id] / units.Hz) ** 0.5 # from elog:1566 and https://en.wikipedia.org/wiki/Johnson%E2%80%93Nyquist_noise (last Eq. in "noise voltage and power" section - logger.status(f'station {station_id} channel {channel_id} noise temperature = {noise_temp_channel}, bandwidth = {self._bandwidth_per_channel[station_id][channel_id]/ units.MHz:.2f} MHz -> Vrms = {self._Vrms_per_channel[station_id][channel_id]/ units.V / units.micro:.2f} muV') + # Calculation of Vrms. For details see from elog:1566 and https://en.wikipedia.org/wiki/Johnson%E2%80%93Nyquist_noise + # (last two Eqs. in "noise voltage and power" section) or our wiki https://nu-radio.github.io/NuRadioMC/NuRadioMC/pages/HDF5_structure.html + + # Bandwidth, i.e., \Delta f in equation + integrated_channel_response = self._integrated_channel_response[station_id][channel_id] + max_amplification = self._max_amplification_per_channel[station_id][channel_id] + + self._Vrms_per_channel[station_id][channel_id] = (noise_temp_channel * 50 * constants.k * integrated_channel_response / units.Hz) ** 0.5 + self._Vrms_efield_per_channel[station_id][channel_id] = self._Vrms_per_channel[station_id][channel_id] / max_amplification / units.m # VEL = 1m + + # for logging + mean_integrated_response = self._integrated_channel_response_normalization[self._station_id][channel_id] + + logger.status(f'Station.channel {station_id}.{channel_id:02d}: noise temperature = {noise_temp_channel} K, ' + f'est. bandwidth = {integrated_channel_response / mean_integrated_response / units.MHz:.2f} MHz, ' + f'max. filter amplification = {max_amplification:.2e} -> Vrms = ' + f'integrated response = {integrated_channel_response / units.MHz:.2e}MHz -> Vrms = ' + f'{self._Vrms_per_channel[station_id][channel_id] / units.mV:.2f} mV -> efield Vrms = {self._Vrms_efield_per_channel[station_id][channel_id] / units.V / units.m / units.micro:.2f}muV/m (assuming VEL = 1m) ') + self._Vrms = next(iter(next(iter(self._Vrms_per_channel.values())).values())) - logger.status('(if same bandwidth for all stations/channels is assumed:) noise temperature = {}, bandwidth = {:.2f} MHz -> Vrms = {:.2f} muV'.format(self._noise_temp, self._bandwidth / units.MHz, self._Vrms / units.V / units.micro)) + elif Vrms is not None: self._Vrms = float(Vrms) * units.V self._noise_temp = None + logger.status(f"Use a fix noise Vrms of {self._Vrms / units.mV:.2f} mV in each channel.") + + for station_id in self._integrated_channel_response: + + for channel_id in self._integrated_channel_response[station_id]: + max_amplification = self._max_amplification_per_channel[station_id][channel_id] + self._Vrms_per_channel[station_id][channel_id] = self._Vrms # to be stored in the hdf5 file + self._Vrms_efield_per_channel[station_id][channel_id] = self._Vrms / max_amplification / units.m # VEL = 1m + + # for logging + integrated_channel_response = self._integrated_channel_response[station_id][channel_id] + mean_integrated_response = self._integrated_channel_response_normalization[self._station_id][channel_id] + + logger.status(f'Station.channel {station_id}.{channel_id:02d}: ' + f'est. bandwidth = {integrated_channel_response / mean_integrated_response / units.MHz:.2f} MHz, ' + f'max. filter amplification = {max_amplification:.2e} ' + f'integrated response = {integrated_channel_response / units.MHz:.2e}MHz ->' + f'efield Vrms = {self._Vrms_efield_per_channel[station_id][channel_id] / units.V / units.m / units.micro:.2f}muV/m (assuming VEL = 1m) ') + else: raise AttributeError(f"noise temperature and Vrms are both set to None") - self._Vrms_efield_per_channel = {} - for station_id in self._bandwidth_per_channel: - self._Vrms_efield_per_channel[station_id] = {} - for channel_id in self._bandwidth_per_channel[station_id]: - self._Vrms_efield_per_channel[station_id][channel_id] = self._Vrms_per_channel[station_id][channel_id] / self._amplification_per_channel[station_id][channel_id] / units.m self._Vrms_efield = next(iter(next(iter(self._Vrms_efield_per_channel.values())).values())) - tmp_cut = float(self._cfg['speedup']['min_efield_amplitude']) - logger.status(f"final Vrms {self._Vrms/units.V:.2g}V corresponds to an efield of {self._Vrms_efield/units.V/units.m/units.micro:.2g} muV/m for a VEL = 1m (amplification factor of system is {amplification:.1f}).\n -> all signals with less then {tmp_cut:.1f} x Vrms_efield = {tmp_cut * self._Vrms_efield/units.m/units.V/units.micro:.2g}muV/m will be skipped") + speed_cut = float(self._cfg['speedup']['min_efield_amplitude']) + logger.status(f"All signals with less then {speed_cut:.1f} x Vrms_efield will be skipped.") self._distance_cut_polynomial = None if self._cfg['speedup']['distance_cut']: @@ -852,7 +903,7 @@ def run(self): self._station = NuRadioReco.framework.station.Station(self._station_id) self._station.set_sim_station(self._sim_station) self._station.get_sim_station().set_station_time(self._evt_time) - + # convert efields to voltages at digitizer if hasattr(self, '_detector_simulation_part1'): # we give the user the opportunity to define a custom detector simulation @@ -955,7 +1006,7 @@ def run(self): channel_ids = self._det.get_channel_ids(self._station.get_id()) Vrms = {} for channel_id in channel_ids: - norm = self._bandwidth_per_channel[self._station.get_id()][channel_id] + norm = self._integrated_channel_response[self._station.get_id()][channel_id] Vrms[channel_id] = self._Vrms_per_channel[self._station.get_id()][channel_id] / (norm / max_freq) ** 0.5 # normalize noise level to the bandwidth its generated for channelGenericNoiseAdder.run(self._evt, self._station, self._det, amplitude=Vrms, min_freq=0 * units.MHz, max_freq=max_freq, type='rayleigh', excluded_channels=self._noiseless_channels[station_id]) @@ -1452,7 +1503,7 @@ def _write_output_file(self, empty=False): positions[channel_id] = self._det.get_relative_position(station_id, channel_id) + self._det.get_absolute_position(station_id) fout["station_{:d}".format(station_id)].attrs['antenna_positions'] = positions 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["station_{:d}".format(station_id)].attrs['bandwidth'] = list(self._integrated_channel_response[station_id].values()) fout.attrs.create("Tnoise", self._noise_temp, dtype=float) fout.attrs.create("Vrms", self._Vrms, dtype=float) diff --git a/documentation/source/NuRadioMC/pages/HDF5_structure.rst b/documentation/source/NuRadioMC/pages/HDF5_structure.rst index 20abb7a84..40fdd4cfb 100644 --- a/documentation/source/NuRadioMC/pages/HDF5_structure.rst +++ b/documentation/source/NuRadioMC/pages/HDF5_structure.rst @@ -45,7 +45,7 @@ The HDF5 files can be thought of as a structured dictionary: - The top level :ref:`attributes `, which can be accessed through ``f.attrs``, contain some top-level information about the simulation. - The :ref:`individual keys ` contain some properties (energy, vertex, ...) for each stored event or shower. -- Finally, the ``station_`` key contains slightly more detailed information (triggers, propagation times, amplitudes...) at the level of individual channels :ref:`for each station `. +- Finally, the ``station_`` key contains slightly more detailed information (triggers, propagation times, amplitudes...) at the level of individual channels :ref:`for each station `. Each station group has its own attributes (``f[station_].attrs``) HDF5 file attributes ____________________ @@ -72,14 +72,26 @@ The top-level attributes can be accessed using ``f.attrs``. These contain: ``start_event_id`` | ``event_id`` of the first event in the file ``trigger_names`` | List of the names of the different triggers simulated ``Tnoise`` | (explicit) noise temperature used in simulation - ``Vrms`` | RMS of the voltage used as thermal noise floor. Determine from ``Tnoise`` and ``bandwidth`` - ``bandwidth`` | Bandwidth of the antennas/detector (for triggering) ``n_samples`` | Samples of the to-be generated antenna signals ``config`` | The (yaml-style) config file used for the simulation ``deposited`` | ``detector`` | The (json-format) detector description used for the simulation ``dt`` | The time resolution, i.e. the inverse of the sampling rate used for the simulation. This is not necessarily the same as the sampling rate of the simulated channels! + +The station-level attributes can be accessed using ``f[station_].attrs``. The first two attributes ``Vrms`` and ``bandwidth`` also exist on the top-level and refer to the corresponding to the first station/channel pair. + + .. _hdf5-station-attrs-table: + + .. csv-table:: HDF5 station attributes + :header: "Key", "Description" + :widths: auto + :delim: | + + ``Vrms`` | RMS of the voltage used as thermal noise floor :math:`v_{n} = (k_{B} \, R \, T \, \Delta f) ^ {0.5}`. See the relevant section "Noise voltage and power" in this `wiki article `_ (last two equations). Determine from ``Tnoise`` and ``bandwidth`` (see below). + ``bandwidth`` | Bandwidth is above equation. Calculated as the integral over the simulated filter response (`filt`) squared: :math:`\Delta f = np.trapz(np.abs(filt) ** 2, ff)`. + ``antenna_positions`` | Relative position of all simulated antennas (channels) + HDF5 file contents __________________ The HDF5 file contains the following items. Listed are the ``key`` and the ``shape`` of each HDF5 dataset, where ``n_events`` is the number of events stored in the file and ``n_showers`` @@ -114,12 +126,12 @@ is the number of showers (which may be larger than the number of events), and `` Station data ____________ In addition, the HDF5 file contains a key for each station in the simulation. -The station contains more detailed information for each station. Some parameters are per event and -some parameters are per shower. See https://doi.org/10.22323/1.395.1231 for a description of how showers relate to events. +The station contains more detailed information for each station. Some parameters are per event and +some parameters are per shower. See https://doi.org/10.22323/1.395.1231 for a description of how showers relate to events. ``m_events`` and ``m_showers`` refer to the number of events and showers that triggered the station. NOTE: The simple table structure of hdf5 files can not capture the complex relation between events and showers in all cases. Some fields can be ambiguous (e.g. `trigger_times` that only lists the last trigger that a shower generated). -For more advanced analyses, please use the ``*.nur`` files. +For more advanced analyses, please use the ``*.nur`` files. The ``event_group_id`` is the same as in the global dictionary. Therefore you can check for one event with an ``event_group_id`` which stations contain the same ``event_group_id`` and retrieve the information, which station triggered, with which amplitude, etc. The same approach works for ``shower_id``.