import numpy as np
import logging
import glob
from NuRadioReco.modules.base.module import register_run
import NuRadioReco.modules.io.NuRadioRecoio
from NuRadioReco.utilities import units
[docs]class channelMeasuredNoiseAdder:
"""
Module that adds measured noise to channel traces
It does so by reading in a set of .nur files, randomly selecting a forced trigger event
and adding the noise from it to the channel waveforms.
The waveforms from the channels in the noise files need to be at least as long as the
waveforms to which the noise is added, so it is recommended to cut them to the right size
first, for example using the channelLengthAdjuster.
"""
def __init__(self):
self.__filenames = None
self.__io = None
self.__random_state = None
self.__max_iterations = None
self.logger = logging.getLogger('NuRadioReco.channelMeasuredNoiseAdder')
self.__noise_data = None
[docs] def begin(self, filenames=None, folder=None, file_pattern="*",
random_seed=None, max_iterations=100, debug=False,
draw_noise_statistics=False, channel_mapping=None, log_level=logging.NOTSET,
restrict_station_id=True, station_id=None, allow_noise_resampling=False,
baseline_substraction=True, allowed_triggers=["FORCE"]):
"""
Set up module parameters
Parameters
----------
filenames: list of strings
List of .nur files containing the measured noise. If None, look for .nur files in "folder".
(Default: None)
folder: str
Only used when "filenames" is None. Directory to search for .nur files matching the "file_pattern"
including subdirectories. (Default: None)
file_pattern: str
Use ``glob.glob(f"{folder}/**/{file_pattern}.nur", recursive=True)`` to search for files. (Default: "*")
random_seed: int, default: None
Seed for the random number generator. By default, no seed is set.
max_iterations: int, default: 100
The module will pick a random event from the noise files, until a suitable event is found
or until the number of iterations exceeds max_iterations. In that case, an error is thrown.
debug: bool, default: False
Set True to get debug output
draw_noise_statistics: boolean, default: False
If true, the values of all samples is stored and a histogram with noise statistics is drawn
be the end() method
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
log_level: logging log level
the log level, default logging.NOTSET (adhere to global logging level)
baseline_substraction: boolean
Option to subtract mean from trace. Set mean_opt=False to remove mean subtraction from trace.
restrict_station_id: bool
Require the station in the noise event to be the same as in the simulated event. (Default: True)
station_id: int
If restrict_station_id is False specify the station id to be used for the noise. If None take first station
in the noise event. (Default: None)
allow_noise_resampling: bool
Allow resampling the noise trace to match the simulated trace. (Default: False)
allowed_triggers: list(str)
List of trigger names which should be used, events with other triggers are not used. (Default: ["FORCE"])
"""
self.logger.setLevel(log_level)
self.__filenames = filenames
if self.__filenames is None:
if folder is None:
err = "Both, \"filenames\" and \"folder\" are None, you have to specify at least one ..."
self.logger.error(err)
raise ValueError(err)
self.__filenames = glob.glob(f"{folder}/**/{file_pattern}.nur", recursive=True)
self.logger.info(f"Found {len(self.__filenames)} noise file(s) ...")
self.__io = NuRadioReco.modules.io.NuRadioRecoio.NuRadioRecoio(self.__filenames)
self.__random_state = np.random.Generator(np.random.Philox(random_seed))
self.__max_iterations = max_iterations
self.__channel_mapping = channel_mapping
self.__baseline_substraction = baseline_substraction
self.__restrict_station_id = restrict_station_id
self.__noise_station_id = station_id
self.__allow_noise_resampling = allow_noise_resampling
self._allowed_triggers = allowed_triggers
if debug:
self.logger.setLevel(logging.DEBUG)
self.logger.debug('Reading noise from {} files containing {} events'.format(len(filenames), self.__io.get_n_events()))
if draw_noise_statistics:
self.__noise_data = []
def __get_noise_channel(self, channel_id):
if self.__channel_mapping is None:
return channel_id
else:
return self.__channel_mapping[channel_id]
[docs] @register_run()
def run(self, event, station, det):
"""
Add measured noise to station channels
Parameters
----------
event: event object
station: station object
det: detector description
"""
noise_station = None
for _ in range(self.__max_iterations):
noise_station = self.get_noise_station(station)
# Get random station from noise file. If we got a suitable station, we continue,
# otherwise we try again
if noise_station is not None:
break
# To avoid infinite loops, if no suitable noise station was found after a number of iterations we raise an error
if noise_station is None:
raise ValueError('Could not find suitable noise event in noise files after {} iterations'.format(self.__max_iterations))
for channel in station.iter_channels():
noise_channel = noise_station.get_channel(self.__get_noise_channel(channel.get_id()))
channel_trace = channel.get_trace()
# if resampling is not desired no channel with wrong sampling rate is selected in get_noise_station()
resampled = False
if noise_channel.get_sampling_rate() != channel.get_sampling_rate():
noise_channel.resample(channel.get_sampling_rate())
resampled = True
noise_trace = noise_channel.get_trace()
if self.__baseline_substraction:
mean = noise_trace.mean()
std = noise_trace.std()
if mean > 0.05 * std:
self.logger.warning((
"The noise trace has an offset/baseline of {:.3f}mV which is more than 5% of the STD of {:.3f}mV. "
"The module corrects for the offset but it might points to an error in the FPN subtraction.").format(mean / units.mV, std / units.mV))
noise_trace -= mean
if len(channel_trace) > len(noise_trace):
err = "{} is shorter than simulated trace. Stop".format("Resampled noise trace" if resampled else "Noise trace")
self.logger.error(err)
raise ValueError(err)
elif len(channel_trace) < len(noise_trace):
self.logger.warn("Noise trace has more samples than the simulated one, clip noise trace at the end ...")
noise_trace = noise_trace[:channel.get_number_of_samples()] # if to long, clip the end
else:
pass
channel_trace += noise_trace
# if draw_noise_statistics == True
if self.__noise_data is not None:
self.__noise_data.append(noise_trace)
[docs] def get_noise_station(self, station):
"""
Returns a random station from the noise files that can be used as a noise sample.
The function selects a random event from the noise files and checks if it is suitable.
If it is, the station is returned, otherwise None is returned. The event is suitable if it
fulfills these criteria:
* It contains a station with the same station ID as the one to which the noise shall be added
* The station does not have a trigger that has triggered.
* The every channel in the station to which the noise shall be added is also present in the station
Parameters
----------
station: Station class
The station to which the noise shall be added
"""
event_i = self.__random_state.integers(0, self.__io.get_n_events())
noise_event = self.__io.get_event_i(event_i)
if self.__restrict_station_id and station.get_id() not in noise_event.get_station_ids():
self.logger.debug('Event {}: No station with ID {} found.'.format(event_i, station.get_id()))
return None
if self.__restrict_station_id:
noise_station = noise_event.get_station(station.get_id())
else:
# If __noise_station_id == None get first station stored in event
noise_station = noise_event.get_station(self.__noise_station_id)
for trigger_name in noise_station.get_triggers():
if trigger_name in self._allowed_triggers:
continue
trigger = noise_station.get_trigger(trigger_name)
if trigger.has_triggered():
self.logger.debug(f'Noise station has triggered ({trigger_name}), reject noise event.')
return None
for channel_id in station.get_channel_ids():
noise_channel_id = self.__get_noise_channel(channel_id)
if noise_channel_id not in noise_station.get_channel_ids():
if channel_id == noise_channel_id:
self.logger.debug('Event {}: Requested channel {} not found'.format(event_i, noise_channel_id))
return None
noise_channel = noise_station.get_channel(noise_channel_id)
channel = station.get_channel(channel_id)
if channel.get_sampling_rate() != noise_channel.get_sampling_rate() and not self.__allow_noise_resampling:
self.logger.debug('Event {}, Channel {}: Different sampling rates, reject noise event.'
.format(event_i, noise_channel_id))
return None
if (noise_channel.get_number_of_samples() / noise_channel.get_sampling_rate() <
channel.get_number_of_samples() / channel.get_sampling_rate()):
# Just add debug message...
self.logger.debug('Event {}, Channel {}: Different sampling rate / trace lenght ratios'
.format(event_i, noise_channel_id))
return noise_station
[docs] def end(self):
"""
End method. Draws a histogram of the noise statistics and fits a
Gaussian distribution to it.
"""
if self.__noise_data is not None:
import matplotlib.pyplot as plt
plt.close('all')
noise_entries = np.array(self.__noise_data)
noise_bins = np.arange(-150, 150, 5.) * units.mV
noise_entries = noise_entries.flatten()
mean = noise_entries.mean()
sigma = noise_entries.std()
fig1, ax1_1 = plt.subplots()
n, bins, pathes = ax1_1.hist(noise_entries / units.mV, bins=noise_bins / units.mV)
ax1_1.plot(noise_bins / units.mV, np.max(n) * np.exp(-0.5 * (noise_bins - mean) ** 2 / sigma ** 2))
ax1_1.grid()
ax1_1.set_xlabel('sample value [mV]')
ax1_1.set_ylabel('entries')
plt.show()