Skip to content

Commit

Permalink
Merge pull request #594 from nu-radio/hotfix-simulation
Browse files Browse the repository at this point in the history 
Renaming bandwidth variable in simulation
  • Loading branch information
@fschlueter
fschlueter authored Jan 19, 2024
2 parents 8542132 + 779573d commit 82ad9f1
Show file tree
Hide file tree
Showing 2 changed files with 100 additions and 37 deletions.
113 changes: 82 additions & 31 deletions NuRadioMC/simulation/simulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -301,66 +302,116 @@ 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)
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)

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']:
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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])
Expand Down Expand Up @@ -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)
Expand Down
24 changes: 18 additions & 6 deletions documentation/source/NuRadioMC/pages/HDF5_structure.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ The HDF5 files can be thought of as a structured dictionary:

- The top level :ref:`attributes <NuRadioMC/pages/HDF5_structure:HDF5 file attributes>`, which can be accessed through ``f.attrs``, contain some top-level information about the simulation.
- The :ref:`individual keys <NuRadioMC/pages/HDF5_structure:HDF5 file contents>` contain some properties (energy, vertex, ...) for each stored event or shower.
- Finally, the ``station_<station_id>`` key contains slightly more detailed information (triggers, propagation times, amplitudes...) at the level of individual channels :ref:`for each station <NuRadioMC/pages/HDF5_structure:Station data>`.
- Finally, the ``station_<station_id>`` key contains slightly more detailed information (triggers, propagation times, amplitudes...) at the level of individual channels :ref:`for each station <NuRadioMC/pages/HDF5_structure:Station data>`. Each station group has its own attributes (``f[station_<station_id>].attrs``)

HDF5 file attributes
____________________
Expand All @@ -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_<station_id>].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 <https://en.wikipedia.org/wiki/Johnson%E2%80%93Nyquist_noise>`_ (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``
Expand Down Expand Up @@ -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``.
Expand Down

0 comments on commit 82ad9f1

Please sign in to comment.