from __future__ import absolute_import, division, print_function
import numpy as np
from radiotools import helper as hp
from radiotools import coordinatesystems as cstrans
from NuRadioMC.SignalGen import askaryan
from NuRadioMC.SignalGen import emitter
from NuRadioReco.utilities import units
from NuRadioMC.utilities import medium
from NuRadioReco.utilities import fft
from NuRadioMC.utilities.earth_attenuation import get_weight
from NuRadioMC.SignalProp import propagation
from numpy.random import Generator, Philox
import h5py
import time
import six
import copy
import json
from scipy import constants
# import detector simulation modules
import NuRadioReco.modules.io.eventWriter
import NuRadioReco.modules.channelSignalReconstructor
import NuRadioReco.modules.electricFieldResampler
import NuRadioReco.modules.channelGenericNoiseAdder
import NuRadioReco.modules.efieldToVoltageConverterPerEfield
import NuRadioReco.modules.efieldToVoltageConverter
import NuRadioReco.modules.channelAddCableDelay
import NuRadioReco.modules.channelResampler
import NuRadioReco.detector.detector as detector
import NuRadioReco.framework.sim_station
import NuRadioReco.framework.electric_field
import NuRadioReco.framework.particle
import NuRadioReco.framework.event
import NuRadioReco.framework.sim_emitter
from NuRadioReco.detector import antennapattern
from NuRadioReco.utilities import geometryUtilities as geo_utl
from NuRadioReco.framework.parameters import channelParameters as chp
from NuRadioReco.framework.parameters import electricFieldParameters as efp
from NuRadioReco.framework.parameters import showerParameters as shp
# parameters describing simulated Monte Carlo particles
from NuRadioReco.framework.parameters import particleParameters as simp
from NuRadioReco.framework.parameters import emitterParameters as ep
# parameters set in the event generator
from NuRadioReco.framework.parameters import generatorAttributes as genattrs
import datetime
import logging
import warnings
from six import iteritems
import yaml
import os
import collections
from NuRadioMC.utilities.Veff import remove_duplicate_triggers
import logging
from NuRadioReco.utilities.logging import LOGGING_STATUS

logger = logging.getLogger("NuRadioMC.simulation")


def pretty_time_delta(seconds):
    seconds = int(seconds)
    days, seconds = divmod(seconds, 86400)
    hours, seconds = divmod(seconds, 3600)
    minutes, seconds = divmod(seconds, 60)
    if days > 0:
        return '%dd%dh%dm%ds' % (days, hours, minutes, seconds)
    elif hours > 0:
        return '%dh%dm%ds' % (hours, minutes, seconds)
    elif minutes > 0:
        return '%dm%ds' % (minutes, seconds)
    else:
        return '%ds' % (seconds,)


def merge_config(user, default):
    if isinstance(user, dict) and isinstance(default, dict):
        for k, v in iteritems(default):
            if k not in user:
                user[k] = v
            else:
                user[k] = merge_config(user[k], v)
    return user


class simulation:

    def __init__(self, inputfilename,
                 outputfilename,
                 detectorfile=None,
                 det=None,
                 det_kwargs={},
                 outputfilenameNuRadioReco=None,
                 debug=False,
                 evt_time=datetime.datetime(2018, 1, 1),
                 config_file=None,
                 log_level=logging.NOTSET,
                 default_detector_station=None,
                 default_detector_channel=None,
                 file_overwrite=False,
                 write_detector=True,
                 event_list=None,
                 log_level_propagation=logging.NOTSET,
                 ice_model=None,
                 **kwargs):
        """
        initialize the NuRadioMC end-to-end simulation

        Parameters
        ----------
        inputfilename: string, or pair
            the path to the hdf5 file containing the list of neutrino events
            alternatively, the data and attributes dictionary can be passed directly to the method
        outputfilename: string
            specify hdf5 output filename.
        detectorfile: string
            path to the json file containing the detector description
        det: detector object
            Pass a detector class object
        det_kwargs: dict
            Pass arguments to the detector (only used when det == None)
        station_id: int
            the station id for which the simulation is performed. Must match a station
            deself._fined in the detector description
        outputfilenameNuRadioReco: string or None
            outputfilename of NuRadioReco detector sim file, this file contains all
            waveforms of the triggered events
            default: None, i.e., no output file will be written which is useful for
            effective volume calculations
        debug: bool
            True activates debug mode, default False
        evt_time: datetime object
            the time of the events, default 1/1/2018
        config_file: string
            path to config file
        log_level: logging.LEVEL
            the log level
        default_detector_station: int or None
            DEPRECATED: Define reference stations in the detector JSON file instead
        default_detector_channel: int or None
            DEPRECATED: Define reference channels in the detector JSON file instead
        file_overwrite: bool
            True allows overwriting of existing files, default False
        write_detector: bool
            If true, the detector description is written into the .nur files along with the events
            default True
        event_list: None or list of ints
            if provided, only the event listed in this list are being simulated
        log_level_propagation: logging.LEVEL
            the log level of the propagation module
        ice_model: medium object (default None)
            allows to specify a custom ice model. This model is used if the config file specifies the ice model as "custom".
        """
        logger.setLevel(log_level)
        if 'write_mode' in kwargs:
            logger.warning('Parameter write_mode is deprecated. Define the output format in the config file instead.')

        self._log_level = log_level
        self._log_level_ray_propagation = log_level_propagation
        config_file_default = os.path.join(os.path.dirname(__file__), 'config_default.yaml')
        logger.status('reading default config from {}'.format(config_file_default))
        with open(config_file_default, 'r') as ymlfile:
            self._cfg = yaml.load(ymlfile, Loader=yaml.FullLoader)
        if config_file is not None:
            logger.status('reading local config overrides from {}'.format(config_file))
            with open(config_file, 'r') as ymlfile:
                local_config = yaml.load(ymlfile, Loader=yaml.FullLoader)
                new_cfg = merge_config(local_config, self._cfg)
                self._cfg = new_cfg

        if self._cfg['seed'] is None:
            # the config seeting None means a random seed. To have the simulation be reproducable, we generate a new
            # random seed once and save this seed to the config setting. If the simulation is rerun, we can get
            # the same random sequence.
            self._cfg['seed'] = np.random.randint(0, 2 ** 32 - 1)

        self._rnd = Generator(Philox(self._cfg['seed']))

        self._outputfilename = outputfilename
        if os.path.exists(self._outputfilename):
            msg = f"hdf5 output file {self._outputfilename} already exists"
            if file_overwrite == False:
                logger.error(msg)
                raise FileExistsError(msg)
            else:
                logger.warning(msg)
        self._detectorfile = detectorfile
        self._n_reflections = int(self._cfg['propagation']['n_reflections'])
        self._outputfilenameNuRadioReco = outputfilenameNuRadioReco
        self._debug = debug
        self._evt_time = evt_time
        self.__write_detector = write_detector
        logger.status("setting event time to {}".format(evt_time))
        self._event_group_list = event_list

        self._antenna_pattern_provider = antennapattern.AntennaPatternProvider()

        # initialize propagation module
        self._prop = propagation.get_propagation_module(self._cfg['propagation']['module'])

        if self._cfg['propagation']['ice_model'] == "custom":
            if ice_model is None:
                logger.error("ice model is set to 'custom' in config file but no custom ice model is provided.")
                raise AttributeError("ice model is set to 'custom' in config file but no custom ice model is provided.")
            self._ice = ice_model
        else:
            self._ice = medium.get_ice_model(self._cfg['propagation']['ice_model'])

        self._mout = collections.OrderedDict()
        self._mout_groups = collections.OrderedDict()
        self._mout_attrs = collections.OrderedDict()

        # Initialize detector
        if det is None:
            logger.status("Detectorfile {}".format(os.path.abspath(self._detectorfile)))
            kwargs = dict(json_filename=self._detectorfile, default_station=default_detector_station,
                              default_channel=default_detector_channel, antenna_by_depth=False)
            kwargs.update(det_kwargs)
            self._det = detector.Detector(**kwargs)
        else:
            self._det = det

        self._det.update(evt_time)

        self._station_ids = self._det.get_station_ids()
        self._event_ids_counter = {}
        for station_id in self._station_ids:
            self._event_ids_counter[station_id] = -1  # we initialize with -1 becaue we increment the counter before we use it the first time

        # print noise information
        logger.status("running with noise {}".format(bool(self._cfg['noise'])))
        logger.status("setting signal to zero {}".format(bool(self._cfg['signal']['zerosignal'])))
        if bool(self._cfg['propagation']['focusing']):
            logger.status("simulating signal amplification due to focusing of ray paths in the firn.")

        # read sampling rate from config (this sampling rate will be used internally)
        self._dt = 1. / (self._cfg['sampling_rate'] * units.GHz)

        if isinstance(inputfilename, str):
            logger.status(f"reading input from {inputfilename}")
            self._inputfilename = inputfilename
            self._read_input_hdf5()  # we read in the full input file into memory at the beginning to limit io to the beginning and end of the run
        else:
            logger.status("getting input on-the-fly")
            self._inputfilename = "on-the-fly"
            self._fin = inputfilename[0]
            self._fin_attrs = inputfilename[1]
            self._fin_stations = {}
        # store all relevant attributes of the input file in a dictionary
        self._generator_info = {}
        for enum_entry in genattrs:
            if enum_entry.name in self._fin_attrs:
                self._generator_info[enum_entry] = self._fin_attrs[enum_entry.name]

        # check if the input file contains events, if not save empty output file (for book keeping) and terminate simulation
        if len(self._fin['xx']) == 0:
            logger.status(f"input file {self._inputfilename} is empty")
            return

        # Perfom a dummy detector simulation to determine how the signals are filtered.
        # This variable stores the integrated channel response for each channel, i.e.
        # the integral of the squared signal chain response over all frequencies, int S21^2 df.
        # For a system without amplification, it is equivalent to the bandwidth of the system.
        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
        self._Vrms = 1
        for iSt, self._station_id in enumerate(self._station_ids):
            self._shower_index = 0
            self._primary_index = 0
            self._evt = NuRadioReco.framework.event.Event(0, self._primary_index)

            self._sampling_rate_detector = self._det.get_sampling_frequency(self._station_id, 0)
#                 logger.warning('internal sampling rate is {:.3g}GHz, final detector sampling rate is {:.3g}GHz'.format(self.get_sampling_rate(), self._sampling_rate_detector))
            self._n_samples = self._det.get_number_of_samples(self._station_id, 0) / self._sampling_rate_detector / self._dt
            self._n_samples = int(np.ceil(self._n_samples / 2.) * 2)  # round to nearest even integer
            self._ff = np.fft.rfftfreq(self._n_samples, self._dt)
            self._tt = np.arange(0, self._n_samples * self._dt, self._dt)

            self._create_sim_station()
            for channel_id in range(self._det.get_number_of_channels(self._station_id)):
                electric_field = NuRadioReco.framework.electric_field.ElectricField([channel_id], self._det.get_relative_position(self._sim_station.get_id(), channel_id))
                trace = np.zeros_like(self._tt)
                trace[self._n_samples // 2] = 100 * units.V  # set a signal that will satisfy any high/low trigger
                trace[self._n_samples // 2 + 1] = -100 * units.V
                electric_field.set_trace(np.array([np.zeros_like(self._tt), trace, trace]), 1. / self._dt)
                electric_field.set_trace_start_time(0)
                electric_field[efp.azimuth] = 0
                electric_field[efp.zenith] = 100 * units.deg
                electric_field[efp.ray_path_type] = 0
                self._sim_station.add_electric_field(electric_field)

            self._station = NuRadioReco.framework.station.Station(self._station_id)
            self._station.set_sim_station(self._sim_station)
            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._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._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._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).")

        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)
                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)

                    # 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()))

        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 = next(iter(next(iter(self._Vrms_efield_per_channel.values())).values()))
        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']:
            coef = self._cfg['speedup']['distance_cut_coefficients']
            self.__distance_cut_polynomial = np.polynomial.polynomial.Polynomial(coef)

            def get_distance_cut(shower_energy):
                if shower_energy <= 0:
                    return 100 * units.m
                return max(100 * units.m, 10 ** self.__distance_cut_polynomial(np.log10(shower_energy)))

            self._get_distance_cut = get_distance_cut

    def run(self):
        """
        run the NuRadioMC simulation
        """
        if len(self._fin['xx']) == 0:
            logger.status(f"writing empty hdf5 output file")
            self._write_output_file(empty=True)
            logger.status(f"terminating simulation")
            return 0
        logger.status(f"Starting NuRadioMC simulation")
        t_start = time.time()
        t_last_update = t_start

        self._channelSignalReconstructor = NuRadioReco.modules.channelSignalReconstructor.channelSignalReconstructor()
        self._eventWriter = NuRadioReco.modules.io.eventWriter.eventWriter()
        efieldToVoltageConverterPerEfield = NuRadioReco.modules.efieldToVoltageConverterPerEfield.efieldToVoltageConverterPerEfield()
        efieldToVoltageConverter = NuRadioReco.modules.efieldToVoltageConverter.efieldToVoltageConverter()
        efieldToVoltageConverter.begin(time_resolution=self._cfg['speedup']['time_res_efieldconverter'])
        channelAddCableDelay = NuRadioReco.modules.channelAddCableDelay.channelAddCableDelay()
        channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder()
        channelGenericNoiseAdder.begin(seed=self._cfg['seed'])
        channelResampler = NuRadioReco.modules.channelResampler.channelResampler()
        electricFieldResampler = NuRadioReco.modules.electricFieldResampler.electricFieldResampler()
        if self._outputfilenameNuRadioReco is not None:
            self._eventWriter.begin(self._outputfilenameNuRadioReco, log_level=self._log_level)
        unique_event_group_ids = np.unique(self._fin['event_group_ids'])
        self._n_showers = len(self._fin['event_group_ids'])
        self._shower_ids = np.array(self._fin['shower_ids'])
        self._shower_index_array = {}  # this array allows to convert the shower id to an index that starts from 0 to be used to access the arrays in the hdf5 file.

        self._raytracer = self._prop(
            self._ice, self._cfg['propagation']['attenuation_model'],
            log_level=self._log_level_ray_propagation,
            n_frequencies_integration=int(self._cfg['propagation']['n_freq']),
            n_reflections=self._n_reflections,
            config=self._cfg,
            detector=self._det
        )
        for shower_index, shower_id in enumerate(self._shower_ids):
            self._shower_index_array[shower_id] = shower_index

        self._create_meta_output_datastructures()

        # check if the same detector was simulated before (then we can save the ray tracing part)
        pre_simulated = self._check_if_was_pre_simulated()

        # Check if vertex_times exists:
        self._check_vertex_times()

        input_time = 0.0
        askaryan_time = 0.0
        rayTracingTime = 0.0
        detSimTime = 0.0
        outputTime = 0.0
        weightTime = 0.0
        distance_cut_time = 0.0

        n_shower_station = len(self._station_ids) * self._n_showers
        iCounter = 0

        # calculate bary centers of station
        self._station_barycenter = np.zeros((len(self._station_ids), 3))
        for iSt, station_id in enumerate(self._station_ids):
            pos = []
            for channel_id in range(self._det.get_number_of_channels(station_id)):
                pos.append(self._det.get_relative_position(station_id, channel_id))
            self._station_barycenter[iSt] = np.mean(np.array(pos), axis=0) + self._det.get_absolute_position(station_id)

        # loop over event groups
        for i_event_group_id, event_group_id in enumerate(unique_event_group_ids):
            logger.debug(f"simulating event group id {event_group_id}")
            if self._event_group_list is not None and event_group_id not in self._event_group_list:
                logger.debug(f"skipping event group {event_group_id} because it is not in the event group list provided to the __init__ function")
                continue
            event_indices = np.atleast_1d(np.squeeze(np.argwhere(self._fin['event_group_ids'] == event_group_id)))

            # the weight calculation is independent of the station, so we do this calculation only once
            # the weight also depends just on the "mother" particle, i.e. the incident neutrino which determines
            # the propability of arriving at our simulation volume. All subsequent showers have the same weight. So
            # we calculate it just once and save it to all subshowers.
            t1 = time.time()

            self._primary_index = event_indices[0]
            # determine if a particle (neutrinos, or a secondary interaction of a neutrino, or surfaec muons) is simulated
            particle_mode = "simulation_mode" not in self._fin_attrs or self._fin_attrs['simulation_mode'] != "emitter"
            self._mout['weights'][event_indices] = np.ones(len(event_indices))  # for a pulser simulation, every event has the same weight
            if particle_mode:
                self._read_input_particle_properties(self._primary_index)  # this sets the self.input_particle for self._primary_index
                # calculate the weight for the primary particle
                self.primary = self.input_particle
                if self._cfg['weights']['weight_mode'] == "existing":
                    if "weights" in self._fin:
                        self._mout['weights'] = self._fin["weights"]
                    else:
                        logger.error("config file specifies to use weights from the input hdf5 file but the input file does not contain this information.")
                elif self._cfg['weights']['weight_mode'] is None:
                    self.primary[simp.weight] = 1.
                else:
                    self.primary[simp.weight] = get_weight(self.primary[simp.zenith],
                                                           self.primary[simp.energy],
                                                           self.primary[simp.flavor],
                                                           mode=self._cfg['weights']['weight_mode'],
                                                           cross_section_type=self._cfg['weights']['cross_section_type'],
                                                           vertex_position=self.primary[simp.vertex],
                                                           phi_nu=self.primary[simp.azimuth])
                # all entries for the event for this primary get the calculated primary's weight
                self._mout['weights'][event_indices] = self.primary[simp.weight]

            weightTime += time.time() - t1
            # skip all events where neutrino weights is zero, i.e., do not
            # simulate neutrino that propagate through the Earth
            if self._mout['weights'][self._primary_index] < self._cfg['speedup']['minimum_weight_cut']:
                logger.debug("neutrino weight is smaller than {}, skipping event".format(self._cfg['speedup']['minimum_weight_cut']))
                continue

            # these quantities get computed to apply the distance cut as a function of shower energies
            # the shower energies of closeby showers will be added as they can constructively interfere
            if self._cfg['speedup']['distance_cut']:
                t_tmp = time.time()
                shower_energies = np.array(self._fin['shower_energies'])[event_indices]
                vertex_positions = np.array([np.array(self._fin['xx'])[event_indices],
                                             np.array(self._fin['yy'])[event_indices],
                                             np.array(self._fin['zz'])[event_indices]]).T
                vertex_distances = np.linalg.norm(vertex_positions - vertex_positions[0], axis=1)
                distance_cut_time += time.time() - t_tmp

            triggered_showers = {}  # this variable tracks which showers triggered a particular station

            # loop over all stations (each station is treated independently)
            for iSt, self._station_id in enumerate(self._station_ids):
                t1 = time.time()
                triggered_showers[self._station_id] = []
                logger.debug(f"simulating station {self._station_id}")

                if self._cfg['speedup']['distance_cut']:
                    # perform a quick cut to reject event group completely if no shower is close enough to the station
                    t_tmp = time.time()
                    vertex_distances_to_station = np.linalg.norm(vertex_positions - self._station_barycenter[iSt], axis=1)
                    distance_cut = self._get_distance_cut(np.sum(shower_energies)) + 100 * units.m  # 100m safety margin is added to account for extent of station around bary center.
                    if vertex_distances_to_station.min() > distance_cut:
                        logger.debug(f"skipping station {self._station_id} because minimal distance {vertex_distances_to_station.min()/units.km:.1f}km > {distance_cut/units.km:.1f}km (shower energy = {shower_energies.max():.2g}eV) bary center of station {self._station_barycenter[iSt]}")
                        distance_cut_time += time.time() - t_tmp
                        iCounter += len(shower_energies)
                        continue
                    distance_cut_time += time.time() - t_tmp

                candidate_station = False
                self._sampling_rate_detector = self._det.get_sampling_frequency(self._station_id, 0)
#                 logger.warning('internal sampling rate is {:.3g}GHz, final detector sampling rate is {:.3g}GHz'.format(self.get_sampling_rate(), self._sampling_rate_detector))
                self._n_samples = self._det.get_number_of_samples(self._station_id, 0) / self._sampling_rate_detector / self._dt
                self._n_samples = int(np.ceil(self._n_samples / 2.) * 2)  # round to nearest even integer
                self._ff = np.fft.rfftfreq(self._n_samples, self._dt)
                self._tt = np.arange(0, self._n_samples * self._dt, self._dt)

                ray_tracing_performed = False
                if 'station_{:d}'.format(self._station_id) in self._fin_stations:
                    ray_tracing_performed = (self._raytracer.get_output_parameters()[0]['name'] in self._fin_stations['station_{:d}'.format(self._station_id)]) and self._was_pre_simulated
                self._evt_tmp = NuRadioReco.framework.event.Event(0, 0)

                if particle_mode:
                    # add the primary particle to the temporary event
                    self._evt_tmp.add_particle(self.primary)

                self._create_sim_station()
                # loop over all showers in event group
                # create output data structure for this channel
                sg = self._create_station_output_structure(len(event_indices), self._det.get_number_of_channels(self._station_id))
                for iSh, self._shower_index in enumerate(event_indices):
                    sg['shower_id'][iSh] = self._shower_ids[self._shower_index]
                    iCounter += 1
#                     if(iCounter % max(1, int(n_shower_station / 100.)) == 0):
                    if (time.time() - t_last_update) > 60 :
                        t_last_update = time.time()
                        eta = pretty_time_delta((time.time() - t_start) * (n_shower_station - iCounter) / iCounter)
                        total_time_sum = input_time + rayTracingTime + detSimTime + outputTime + weightTime + distance_cut_time  # askaryan time is part of the ray tracing time, so it is not counted here.
                        total_time = time.time() - t_start
                        tmp_att = 0
                        if total_time > 0:
                            logger.status(
                                "processing event group {}/{} and shower {}/{} ({} showers triggered) = {:.1f}%, ETA {}, time consumption: ray tracing = {:.0f}%, askaryan = {:.0f}%, detector simulation = {:.0f}% reading input = {:.0f}%, calculating weights = {:.0f}%, distance cut {:.0f}%, unaccounted = {:.0f}% ".format(
                                    i_event_group_id,
                                    len(unique_event_group_ids),
                                    iCounter,
                                    n_shower_station,
                                    np.sum(self._mout['triggered']),
                                    100. * iCounter / n_shower_station,
                                    eta,
                                    100. * (rayTracingTime - askaryan_time) / total_time,
                                    100. * askaryan_time / total_time,
                                    100. * detSimTime / total_time,
                                    100.*input_time / total_time,
                                    100. * weightTime / total_time,
                                    100 * distance_cut_time / total_time,
                                    100 * (total_time - total_time_sum) / total_time))

                    self._read_input_shower_properties()
                    if particle_mode:
                        logger.debug(f"simulating shower {self._shower_index}: {self._fin['shower_type'][self._shower_index]} with E = {self._fin['shower_energies'][self._shower_index]/units.eV:.2g}eV")
                    x1 = self._shower_vertex  # the interaction point

                    if self._cfg['speedup']['distance_cut']:
                        t_tmp = time.time()
                        # calculate the sum of shower energies for all showers within self._cfg['speedup']['distance_cut_sum_length']
                        mask_shower_sum = np.abs(vertex_distances - vertex_distances[iSh]) < self._cfg['speedup']['distance_cut_sum_length']
                        shower_energy_sum = np.sum(shower_energies[mask_shower_sum])
                        # quick speedup cut using barycenter of station as position
                        distance_to_station = np.linalg.norm(x1 - self._station_barycenter[iSt])
                        distance_cut = self._get_distance_cut(shower_energy_sum) + 100 * units.m  # 100m safety margin is added to account for extent of station around bary center.
                        logger.debug(f"calculating distance cut. Current event has energy {self._fin['shower_energies'][self._shower_index]:.4g}, it is event number {iSh} and {np.sum(mask_shower_sum)} are within {self._cfg['speedup']['distance_cut_sum_length']/units.m:.1f}m -> {shower_energy_sum:.4g}")
                        if distance_to_station > distance_cut:
                            logger.debug(f"skipping station {self._station_id} because distance {distance_to_station/units.km:.1f}km > {distance_cut/units.km:.1f}km (shower energy = {self._fin['shower_energies'][self._shower_index]:.2g}eV) between vertex {x1} and bary center of station {self._station_barycenter[iSt]}")
                            distance_cut_time += time.time() - t_tmp
                            continue
                        distance_cut_time += time.time() - t_tmp

                    # skip vertices not in fiducial volume. This is required because 'mother' events are added to the event list
                    # if daugthers (e.g. tau decay) have their vertex in the fiducial volume
                    if not self._is_in_fiducial_volume(self._shower_vertex):
                        logger.debug(f"event is not in fiducial volume, skipping simulation {self._fin['xx'][self._shower_index]}, "
                                     f"{self._fin['yy'][self._shower_index]}, {self._fin['zz'][self._shower_index]}")
                        continue

                    # for special cases where only EM or HAD showers are simulated, skip all events that don't fulfill this criterion
                    if self._cfg['signal']['shower_type'] == "em":
                        if self._fin['shower_type'][self._shower_index] != "em":
                            continue
                    if self._cfg['signal']['shower_type'] == "had":
                        if self._fin['shower_type'][self._shower_index] != "had":
                            continue

                    if particle_mode:
                        self._create_sim_shower()  # create sim shower
                        self._evt_tmp.add_sim_shower(self._sim_shower)
                    else:
                        emitter_obj = NuRadioReco.framework.sim_emitter.SimEmitter(self._shower_index)  # shower_id is equivalent to emitter_id in this case
                        emitter_obj[ep.position] = np.array([self._fin['xx'][self._primary_index], self._fin['yy'][self._primary_index], self._fin['zz'][self._primary_index]])
                        emitter_obj[ep.model] = self._fin['emitter_model'][self._primary_index]
                        emitter_obj[ep.amplitude] = self._fin['emitter_amplitudes'][self._primary_index]
                        for key in ep:
                            if not emitter_obj.has_parameter(key):
                                if 'emitter_' + key.name in self._fin:
                                    emitter_obj[key] = self._fin['emitter_' + key.name][self._primary_index]
                        self._evt_tmp.add_sim_emitter(emitter_obj)

                    # generate unique and increasing event id per station
                    self._event_ids_counter[self._station_id] += 1
                    self._event_id = self._event_ids_counter[self._station_id]

                    # be careful, zenith/azimuth angle always refer to where the neutrino came from,
                    # i.e., opposite to the direction of propagation. We need the propagation direction here,
                    # so we multiply the shower axis with '-1'
                    if 'zeniths' in self._fin:
                        self._shower_axis = -1 * hp.spherical_to_cartesian(self._fin['zeniths'][self._shower_index], self._fin['azimuths'][self._shower_index])
                    else:
                        self._shower_axis = np.array([0, 0, 1])

                    # calculate correct Cherenkov angle for ice density at vertex position
                    n_index = self._ice.get_index_of_refraction(x1)
                    cherenkov_angle = np.arccos(1. / n_index)

                    # first step: perform raytracing to see if solution exists
                    t2 = time.time()
#                     input_time += (time.time() - t1)

                    for channel_id in range(self._det.get_number_of_channels(self._station_id)):
                        x2 = self._det.get_relative_position(self._station_id, channel_id) + self._det.get_absolute_position(self._station_id)
                        logger.debug(f"simulating channel {channel_id} at {x2}")

                        if self._cfg['speedup']['distance_cut']:
                            t_tmp = time.time()
                            distance_cut = self._get_distance_cut(shower_energy_sum)
                            distance = np.linalg.norm(x1 - x2)

                            if distance > distance_cut:
                                logger.debug('A distance speed up cut has been applied')
                                logger.debug('Shower energy: {:.2e} eV'.format(self._fin['shower_energies'][self._shower_index] / units.eV))
                                logger.debug('Distance cut: {:.2f} m'.format(distance_cut / units.m))
                                logger.debug('Distance to vertex: {:.2f} m'.format(distance / units.m))
                                distance_cut_time += time.time() - t_tmp
                                continue
                            distance_cut_time += time.time() - t_tmp

                        self._raytracer.set_start_and_end_point(x1, x2)
                        self._raytracer.use_optional_function('set_shower_axis', self._shower_axis)
                        if pre_simulated and ray_tracing_performed and not self._cfg['speedup']['redo_raytracing']:  # check if raytracing was already performed
                            if self._cfg['propagation']['module'] == 'radiopropa':
                                logger.error('Presimulation can not be used with the radiopropa ray tracer module')
                                raise Exception('Presimulation can not be used with the radiopropa ray tracer module')
                            sg_pre = self._fin_stations["station_{:d}".format(self._station_id)]
                            ray_tracing_solution = {}
                            for output_parameter in self._raytracer.get_output_parameters():
                                ray_tracing_solution[output_parameter['name']] = sg_pre[output_parameter['name']][self._shower_index, channel_id]
                            self._raytracer.set_solution(ray_tracing_solution)
                        else:
                            self._raytracer.find_solutions()

                        if not self._raytracer.has_solution():
                            logger.debug("event {} and station {}, channel {} does not have any ray tracing solution ({} to {})".format(
                                self._event_group_id, self._station_id, channel_id, x1, x2))
                            continue
                        delta_Cs = []
                        viewing_angles = []
                        # loop through all ray tracing solution
                        for iS in range(self._raytracer.get_number_of_solutions()):
                            for key, value in self._raytracer.get_raytracing_output(iS).items():
                                sg[key][iSh, channel_id, iS] = value
                            self._launch_vector = self._raytracer.get_launch_vector(iS)
                            sg['launch_vectors'][iSh, channel_id, iS] = self._launch_vector
                            # calculates angle between shower axis and launch vector
                            viewing_angle = hp.get_angle(self._shower_axis, self._launch_vector)
                            viewing_angles.append(viewing_angle)
                            delta_C = (viewing_angle - cherenkov_angle)
                            logger.debug('solution {} {}: viewing angle {:.1f} = delta_C = {:.1f}'.format(
                                iS, propagation.solution_types[self._raytracer.get_solution_type(iS)], viewing_angle / units.deg, (viewing_angle - cherenkov_angle) / units.deg))
                            delta_Cs.append(delta_C)

                        # discard event if delta_C (angle off cherenkov cone) is too large
                        if min(np.abs(delta_Cs)) > self._cfg['speedup']['delta_C_cut']:
                            logger.debug('delta_C too large, event unlikely to be observed, skipping event')
                            continue

                        n = self._raytracer.get_number_of_solutions()
                        for iS in range(n):  # loop through all ray tracing solution
                            # skip individual channels where the viewing angle difference is too large
                            # discard event if delta_C (angle off cherenkov cone) is too large
                            if np.abs(delta_Cs[iS]) > self._cfg['speedup']['delta_C_cut']:
                                logger.debug('delta_C too large, ray tracing solution unlikely to be observed, skipping event')
                                continue

                            if pre_simulated and ray_tracing_performed and not self._cfg['speedup']['redo_raytracing']:
                                sg_pre = self._fin_stations["station_{:d}".format(self._station_id)]

                                wave_propagation_distance = sg_pre['travel_distances'][self._shower_index, channel_id, iS]
                                wave_propagation_time = sg_pre['travel_times'][self._shower_index, channel_id, iS]
                            else:
                                wave_propagation_distance = self._raytracer.get_path_length(iS)  # calculate path length
                                wave_propagation_time = self._raytracer.get_travel_time(iS)  # calculate travel time
                                if wave_propagation_distance is None or wave_propagation_time is None:
                                    continue

                            sg['travel_distances'][iSh, channel_id, iS] = wave_propagation_distance
                            sg['travel_times'][iSh, channel_id, iS] = wave_propagation_time

                            self._launch_vector = self._raytracer.get_launch_vector(iS)
                            receive_vector = self._raytracer.get_receive_vector(iS)
                            # save receive vector
                            sg['receive_vectors'][iSh, channel_id, iS] = receive_vector
                            zenith, azimuth = hp.cartesian_to_spherical(*receive_vector)

                            # get neutrino pulse from Askaryan module
                            t_ask = time.time()

                            if "simulation_mode" not in self._fin_attrs or self._fin_attrs['simulation_mode'] == "neutrino":
                                # first consider in-ice showers
                                kwargs = {}
                                # if the input file specifies a specific shower realization, use that realization
                                if self._cfg['signal']['model'] in ["ARZ2019", "ARZ2020"] and "shower_realization_ARZ" in self._fin:
                                    kwargs['iN'] = self._fin['shower_realization_ARZ'][self._shower_index]
                                    logger.debug(f"reusing shower {kwargs['iN']} ARZ shower library")
                                elif self._cfg['signal']['model'] == "Alvarez2009" and "shower_realization_Alvarez2009" in self._fin:
                                    kwargs['k_L'] = self._fin['shower_realization_Alvarez2009'][self._shower_index]
                                    logger.debug(f"reusing k_L parameter of Alvarez2009 model of k_L = {kwargs['k_L']:.4g}")
                                else:
                                    # check if the shower was already simulated (e.g. for a different channel or ray tracing solution)
                                    if self._cfg['signal']['model'] in ["ARZ2019", "ARZ2020"]:
                                        if self._sim_shower.has_parameter(shp.charge_excess_profile_id):
                                            kwargs = {'iN': self._sim_shower.get_parameter(shp.charge_excess_profile_id)}
                                    if self._cfg['signal']['model'] == "Alvarez2009":
                                        if self._sim_shower.has_parameter(shp.k_L):
                                            kwargs = {'k_L': self._sim_shower.get_parameter(shp.k_L)}
                                            logger.debug(f"reusing k_L parameter of Alvarez2009 model of k_L = {kwargs['k_L']:.4g}")

                                spectrum, additional_output = askaryan.get_frequency_spectrum(self._fin['shower_energies'][self._shower_index], viewing_angles[iS],
                                                self._n_samples, self._dt, self._fin['shower_type'][self._shower_index], n_index, wave_propagation_distance,
                                                self._cfg['signal']['model'], seed=self._cfg['seed'], full_output=True, **kwargs)
                                # save shower realization to SimShower and hdf5 file
                                if self._cfg['signal']['model'] in ["ARZ2019", "ARZ2020"]:
                                    if 'shower_realization_ARZ' not in self._mout:
                                        self._mout['shower_realization_ARZ'] = np.zeros(self._n_showers)
                                    if not self._sim_shower.has_parameter(shp.charge_excess_profile_id):
                                        self._sim_shower.set_parameter(shp.charge_excess_profile_id, additional_output['iN'])
                                        self._mout['shower_realization_ARZ'][self._shower_index] = additional_output['iN']
                                        logger.debug(f"setting shower profile for ARZ shower library to i = {additional_output['iN']}")
                                if self._cfg['signal']['model'] == "Alvarez2009":
                                    if 'shower_realization_Alvarez2009' not in self._mout:
                                        self._mout['shower_realization_Alvarez2009'] = np.zeros(self._n_showers)
                                    if not self._sim_shower.has_parameter(shp.k_L):
                                        self._sim_shower.set_parameter(shp.k_L, additional_output['k_L'])
                                        self._mout['shower_realization_Alvarez2009'][self._shower_index] = additional_output['k_L']
                                        logger.debug(f"setting k_L parameter of Alvarez2009 model to k_L = {additional_output['k_L']:.4g}")
                                askaryan_time += (time.time() - t_ask)

                                polarization_direction_onsky = self._calculate_polarization_vector()
                                cs_at_antenna = cstrans.cstrafo(*hp.cartesian_to_spherical(*receive_vector))
                                polarization_direction_at_antenna = cs_at_antenna.transform_from_onsky_to_ground(polarization_direction_onsky)
                                logger.debug('receive zenith {:.0f} azimuth {:.0f} polarization on sky {:.2f} {:.2f} {:.2f}, on ground @ antenna {:.2f} {:.2f} {:.2f}'.format(
                                    zenith / units.deg, azimuth / units.deg, polarization_direction_onsky[0],
                                    polarization_direction_onsky[1], polarization_direction_onsky[2],
                                    *polarization_direction_at_antenna))
                                sg['polarization'][iSh, channel_id, iS] = polarization_direction_at_antenna
                                eR, eTheta, ePhi = np.outer(polarization_direction_onsky, spectrum)

                            elif self._fin_attrs['simulation_mode'] == "emitter":
                                # NuRadioMC also supports the simulation of emitters. In this case, the signal model specifies the electric field polarization
                                amplitude = self._fin['emitter_amplitudes'][self._shower_index]
                                emitter_model = self._fin['emitter_model'][self._shower_index]
                                emitter_kwargs = {}
                                emitter_kwargs["launch_vector"] = self._launch_vector
                                for key in self._fin.keys():
                                    if key not in ['emitter_amplitudes', 'emitter_model']:
                                        if key.startswith("emitter_"):
                                            emitter_kwargs[key[8:]] = self._fin[key][self._shower_index]

                                if emitter_model.startswith("efield_"):
                                    if emitter_model == "efield_idl1_spice":
                                        if "emitter_realization" in self._fin:
                                            emitter_kwargs['iN'] = self._fin['emitter_realization'][self._shower_index]
                                        elif emitter_obj.has_parameter(ep.realization_id):
                                            emitter_kwargs['iN'] = emitter_obj.get_parameter(ep.realization_id)
                                        else:
                                            emitter_kwargs['rnd'] = self._rnd

                                    (eR, eTheta, ePhi), additional_output = emitter.get_frequency_spectrum(amplitude, self._n_samples, self._dt, emitter_model, **emitter_kwargs, full_output=True)
                                    if emitter_model == "efield_idl1_spice":
                                        if 'emitter_realization' not in self._mout:
                                            self._mout['emitter_realization'] = np.zeros(self._n_showers)
                                        if not emitter_obj.has_parameter(ep.realization_id):
                                            emitter_obj.set_parameter(ep.realization_id, additional_output['iN'])
                                            self._mout['emitter_realization'][self._shower_index] = additional_output['iN']
                                            logger.debug(f"setting emitter realization to i = {additional_output['iN']}")
                                else:
                                    # the emitter fuction returns the voltage output of the pulser. We need to convole with the antenna response of the emitting antenna
                                    # to obtain the emitted electric field.
                                    # get emitting antenna properties
                                    antenna_model = self._fin['emitter_antenna_type'][self._shower_index]
                                    antenna_pattern = self._antenna_pattern_provider.load_antenna_pattern(antenna_model)
                                    ori = [self._fin['emitter_orientation_theta'][self._shower_index], self._fin['emitter_orientation_phi'][self._shower_index],
                                           self._fin['emitter_rotation_theta'][self._shower_index], self._fin['emitter_rotation_phi'][self._shower_index]]

                                    # source voltage given to the emitter
                                    voltage_spectrum_emitter = emitter.get_frequency_spectrum(amplitude, self._n_samples, self._dt,
                                                                                              emitter_model, **emitter_kwargs)
                                    # convolve voltage output with antenna response to obtain emitted electric field
                                    frequencies = np.fft.rfftfreq(self._n_samples, d=self._dt)
                                    zenith_emitter, azimuth_emitter = hp.cartesian_to_spherical(*self._launch_vector)
                                    VEL = antenna_pattern.get_antenna_response_vectorized(frequencies, zenith_emitter, azimuth_emitter, *ori)
                                    c = constants.c * units.m / units.s
                                    eTheta = VEL['theta'] * (-1j) * voltage_spectrum_emitter * frequencies * n_index / c
                                    ePhi = VEL['phi'] * (-1j) * voltage_spectrum_emitter * frequencies * n_index / c
                                    eR = np.zeros_like(eTheta)
                                # rescale amplitudes by 1/R, for emitters this is not part of the "SignalGen" class
                                eTheta *= 1 / wave_propagation_distance
                                ePhi *= 1 / wave_propagation_distance
                            else:
                                logger.error(f"simulation mode {self._fin_attrs['simulation_mode']} unknown.")
                                raise AttributeError(f"simulation mode {self._fin_attrs['simulation_mode']} unknown.")

                            if self._debug:
                                from matplotlib import pyplot as plt
                                fig, (ax, ax2) = plt.subplots(1, 2)
                                ax.plot(self._ff, np.abs(eTheta) / units.micro / units.V * units.m)
                                ax2.plot(self._tt, fft.freq2time(eTheta, 1. / self._dt) / units.micro / units.V * units.m)
                                ax2.set_ylabel("amplitude [$\mu$V/m]")
                                fig.tight_layout()
                                fig.suptitle("$E_C$ = {:.1g}eV $\Delta \Omega$ = {:.1f}deg, R = {:.0f}m".format(
                                    self._fin['shower_energies'][self._shower_index], viewing_angles[iS], wave_propagation_distance))
                                fig.subplots_adjust(top=0.9)
                                plt.show()

                            electric_field = NuRadioReco.framework.electric_field.ElectricField([channel_id],
                                                position=self._det.get_relative_position(self._sim_station.get_id(), channel_id),
                                                shower_id=self._shower_ids[self._shower_index], ray_tracing_id=iS)
                            if iS is None:
                                a = 1 / 0
                            electric_field.set_frequency_spectrum(np.array([eR, eTheta, ePhi]), 1. / self._dt)
                            electric_field = self._raytracer.apply_propagation_effects(electric_field, iS)
                            # Trace start time is equal to the interaction time relative to the first
                            # interaction plus the wave travel time.
                            if hasattr(self, '_vertex_time'):
                                trace_start_time = self._vertex_time + wave_propagation_time
                            else:
                                trace_start_time = wave_propagation_time

                            # We shift the trace start time so that the trace time matches the propagation time.
                            # The centre of the trace corresponds to the instant when the signal from the shower
                            # vertex arrives at the observer. The next line makes sure that the centre time
                            # of the trace is equal to vertex_time + wave_propagation_time (wave propagation time)
                            trace_start_time -= 0.5 * electric_field.get_number_of_samples() / electric_field.get_sampling_rate()

                            electric_field.set_trace_start_time(trace_start_time)
                            electric_field[efp.azimuth] = azimuth
                            electric_field[efp.zenith] = zenith
                            electric_field[efp.ray_path_type] = propagation.solution_types[self._raytracer.get_solution_type(iS)]
                            electric_field[efp.nu_vertex_distance] = sg['travel_distances'][iSh, channel_id, iS]
                            electric_field[efp.nu_vertex_propagation_time] = sg['travel_times'][iSh, channel_id, iS]
                            electric_field[efp.nu_viewing_angle] = viewing_angles[iS]
                            self._sim_station.add_electric_field(electric_field)

                            # apply a simple threshold cut to speed up the simulation,
                            # application of antenna response will just decrease the
                            # signal amplitude
                            if np.max(np.abs(electric_field.get_trace())) > float(self._cfg['speedup']['min_efield_amplitude']) * self._Vrms_efield_per_channel[self._station_id][channel_id]:
                                candidate_station = True
                        # end of ray tracing solutions loop
                    t3 = time.time()
                    rayTracingTime += t3 - t2
                    # end of channels loop
                # end of showers loop
                # now perform first part of detector simulation -> convert each efield to voltage
                # (i.e. apply antenna response) and apply additional simulation of signal chain (such as cable delays,
                # amp response etc.)
                if not candidate_station:
                    logger.debug("electric field amplitude too small in all channels, skipping to next event")
                    continue
                t1 = time.time()
                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
                    self._detector_simulation_part1()
                else:
                    efieldToVoltageConverterPerEfield.run(self._evt, self._station, self._det)  # convolve efield with antenna pattern
                    self._detector_simulation_filter_amp(self._evt, self._station.get_sim_station(), self._det)
                    channelAddCableDelay.run(self._evt, self._sim_station, self._det)

                if self._cfg['speedup']['amp_per_ray_solution']:
                    self._channelSignalReconstructor.run(self._evt, self._station.get_sim_station(), self._det)
                    for channel in self._station.get_sim_station().iter_channels():
                        tmp_index = np.argwhere(event_indices == self._get_shower_index(channel.get_shower_id()))[0]
                        sg['max_amp_shower_and_ray'][tmp_index, channel.get_id(), channel.get_ray_tracing_solution_id()] = channel.get_parameter(chp.maximum_amplitude_envelope)
                        sg['time_shower_and_ray'][tmp_index, channel.get_id(), channel.get_ray_tracing_solution_id()] = channel.get_parameter(chp.signal_time)
                start_times = []
                channel_identifiers = []
                for channel in self._sim_station.iter_channels():
                    channel_identifiers.append(channel.get_unique_identifier())
                    start_times.append(channel.get_trace_start_time())
                start_times = np.array(start_times)
                start_times_sort = np.argsort(start_times)
                delta_start_times = start_times[start_times_sort][1:] - start_times[start_times_sort][:-1]  # this array is sorted in time
                split_event_time_diff = float(self._cfg['split_event_time_diff'])
                iSplit = np.atleast_1d(np.squeeze(np.argwhere(delta_start_times > split_event_time_diff)))
#                 print(f"start times {start_times}")
#                 print(f"sort array {start_times_sort}")
#                 print(f"delta times {delta_start_times}")
#                 print(f"split at indices {iSplit}")
                n_sub_events = len(iSplit) + 1
                if n_sub_events > 1:
                    logger.info(f"splitting event group id {self._event_group_id} into {n_sub_events} sub events")

                tmp_station = copy.deepcopy(self._station)
                event_group_has_triggered = False
                for iEvent in range(n_sub_events):
                    iStart = 0
                    iStop = len(channel_identifiers)
                    if n_sub_events > 1:
                        if iEvent > 0:
                            iStart = iSplit[iEvent - 1] + 1
                    if iEvent < n_sub_events - 1:
                        iStop = iSplit[iEvent] + 1
                    indices = start_times_sort[iStart: iStop]
                    if n_sub_events > 1:
                        tmp = ""
                        for start_time in start_times[indices]:
                            tmp += f"{start_time/units.ns:.0f}, "
                        tmp = tmp[:-2] + " ns"
                        logger.info(f"creating event {iEvent} of event group {self._event_group_id} ranging rom {iStart} to {iStop} with indices {indices} corresponding to signal times of {tmp}")
                    self._evt = NuRadioReco.framework.event.Event(self._event_group_id, iEvent)  # create new event

                    if particle_mode:
                        # add MC particles that belong to this (sub) event to event structure
                        # add only primary for now, since full interaction chain is not typically in the input hdf5s
                        self._evt.add_particle(self.primary)
                    # copy over generator information from temporary event to event
                    self._evt._generator_info = self._generator_info

                    self._station = NuRadioReco.framework.station.Station(self._station_id)
                    sim_station = NuRadioReco.framework.sim_station.SimStation(self._station_id)
                    sim_station.set_is_neutrino()
                    tmp_sim_station = tmp_station.get_sim_station()
                    self._shower_ids_of_sub_event = []
                    for iCh in indices:
                        ch_uid = channel_identifiers[iCh]
                        shower_id = ch_uid[1]
                        if shower_id not in self._shower_ids_of_sub_event:
                            self._shower_ids_of_sub_event.append(shower_id)
                        sim_station.add_channel(tmp_sim_station.get_channel(ch_uid))
                        efield_uid = ([ch_uid[0]], ch_uid[1], ch_uid[2])  # the efield unique identifier has as first parameter an array of the channels it is valid for
                        for efield in tmp_sim_station.get_electric_fields():
                            if efield.get_unique_identifier() == efield_uid:
                                sim_station.add_electric_field(efield)

                    if particle_mode:
                        # add showers that contribute to this (sub) event to event structure
                        for shower_id in self._shower_ids_of_sub_event:
                            self._evt.add_sim_shower(self._evt_tmp.get_sim_shower(shower_id))
                    else:
                        for shower_id in self._shower_ids_of_sub_event:
                            self._evt.add_sim_emitter(self._evt_tmp.get_sim_emitter(shower_id))
                    self._station.set_sim_station(sim_station)
                    self._station.set_station_time(self._evt_time)
                    self._evt.set_station(self._station)
                    if bool(self._cfg['signal']['zerosignal']):
                        self._increase_signal(None, 0)

                    logger.debug("performing detector simulation")
                    if hasattr(self, '_detector_simulation_part2'):
                        # we give the user the opportunity to specify a custom detector simulation module sequence
                        # which might be needed for certain analyses
                        self._detector_simulation_part2()
                    else:
                        # start detector simulation
                        efieldToVoltageConverter.run(self._evt, self._station, self._det)  # convolve efield with antenna pattern
                        # downsample trace to internal simulation sampling rate (the efieldToVoltageConverter upsamples the trace to
                        # 20 GHz by default to achive a good time resolution when the two signals from the two signal paths are added)
                        channelResampler.run(self._evt, self._station, self._det, sampling_rate=1. / self._dt)

                        if self._is_simulate_noise():
                            max_freq = 0.5 / self._dt
                            channel_ids = self._det.get_channel_ids(self._station.get_id())
                            Vrms = {}
                            for channel_id in channel_ids:
                                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])

                        self._detector_simulation_filter_amp(self._evt, self._station, self._det)

                        self._detector_simulation_trigger(self._evt, self._station, self._det)
                    if not self._station.has_triggered():
                        continue

                    event_group_has_triggered = True
                    triggered_showers[self._station_id].extend(self._get_shower_index(self._shower_ids_of_sub_event))
                    self._calculate_signal_properties()

                    global_shower_indices = self._get_shower_index(self._shower_ids_of_sub_event)
                    local_shower_index = np.atleast_1d(np.squeeze(np.argwhere(np.isin(event_indices, global_shower_indices, assume_unique=True))))
                    self._save_triggers_to_hdf5(sg, local_shower_index, global_shower_indices)
                    if self._outputfilenameNuRadioReco is not None:
                        # downsample traces to detector sampling rate to save file size
                        channelResampler.run(self._evt, self._station, self._det, sampling_rate=self._sampling_rate_detector)
                        channelResampler.run(self._evt, self._station.get_sim_station(), self._det, sampling_rate=self._sampling_rate_detector)
                        electricFieldResampler.run(self._evt, self._station.get_sim_station(), self._det, sampling_rate=self._sampling_rate_detector)

                        output_mode = {'Channels': self._cfg['output']['channel_traces'],
                                       'ElectricFields': self._cfg['output']['electric_field_traces'],
                                       'SimChannels': self._cfg['output']['sim_channel_traces'],
                                       'SimElectricFields': self._cfg['output']['sim_electric_field_traces']}
                        if self.__write_detector:
                            self._eventWriter.run(self._evt, self._det, mode=output_mode)
                        else:
                            self._eventWriter.run(self._evt, mode=output_mode)
                        logger.debug("WRITING EVENT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
                # end sub events loop

                # add local sg array to output data structure if any
                if event_group_has_triggered:
                    if self._station_id not in self._mout_groups:
                        self._mout_groups[self._station_id] = {}
                    for key in sg:
                        if key not in self._mout_groups[self._station_id]:
                            self._mout_groups[self._station_id][key] = list(sg[key])
                        else:
                            self._mout_groups[self._station_id][key].extend(sg[key])

                detSimTime += time.time() - t1

            # end station loop

        # end event group loop

        # Create trigger structures if there are no triggering events.
        # This is done to ensure that files with no triggering n_events
        # merge properly.
#         self._create_empty_multiple_triggers()

        # save simulation run in hdf5 format (only triggered events)
        t5 = time.time()
        self._write_output_file()
        if self._outputfilenameNuRadioReco is not None:
            self._eventWriter.end()
            logger.debug("closing nur file")

        if "simulation_mode" not in self._fin_attrs or self._fin_attrs['simulation_mode'] == "neutrino": # only calcualte Veff for neutrino simulations
            try:
                self.calculate_Veff()
            except:
                logger.error("error in calculating effective volume")

        t_total = time.time() - t_start
        outputTime = time.time() - t5

        output_NuRadioRecoTime = "Timing of NuRadioReco modules \n"
        ts = []
        for iM, (name, instance, kwargs) in enumerate(self._evt.iter_modules(self._station.get_id())):
            ts.append(instance.run.time[instance])
        ttot = np.sum(np.array(ts))
        for i, (name, instance, kwargs) in enumerate(self._evt.iter_modules(self._station.get_id())):
            t = pretty_time_delta(ts[i])
            trel = 100.*ts[i] / ttot
            output_NuRadioRecoTime += f"{name}: {t} {trel:.1f}%\n"
        logger.status(output_NuRadioRecoTime)

        logger.status("{:d} events processed in {} = {:.2f}ms/event ({:.1f}% input, {:.1f}% ray tracing, {:.1f}% askaryan, {:.1f}% detector simulation, {:.1f}% output, {:.1f}% weights calculation)".format(self._n_showers,
                                                                                         pretty_time_delta(t_total), 1.e3 * t_total / self._n_showers,
                                                                                         100 * input_time / t_total,
                                                                                         100 * (rayTracingTime - askaryan_time) / t_total,
                                                                                         100 * askaryan_time / t_total,
                                                                                         100 * detSimTime / t_total,
                                                                                         100 * outputTime / t_total,
                                                                                         100 * weightTime / t_total))
        triggered = remove_duplicate_triggers(self._mout['triggered'], self._fin['event_group_ids'])
        n_triggered = np.sum(triggered)
        return n_triggered

    def _calculate_emitter_output(self):
        pass

    def _get_shower_index(self, shower_id):
        if hasattr(shower_id, "__len__"):
            return np.array([self._shower_index_array[x] for x in shower_id])
        else:
            return self._shower_index_array[shower_id]

    def _is_simulate_noise(self):
        """
        returns True if noise should be added
        """
        return bool(self._cfg['noise'])

    def _is_in_fiducial_volume(self, pos):
        """ Checks if pos is in fiducial volume """

        for check_attr in ['fiducial_zmin', 'fiducial_zmax']:
            if not check_attr in self._fin_attrs:
                logger.warning("Fiducial volume not defined. Return True")
                return True

        pos = copy.deepcopy(pos) - np.array([self._fin_attrs.get("x0", 0), self._fin_attrs.get("y0", 0), 0])

        if not (self._fin_attrs["fiducial_zmin"] < pos[2] < self._fin_attrs["fiducial_zmax"]):
            return False

        if "fiducial_rmax" in self._fin_attrs:
            radius = np.sqrt(pos[0] ** 2 + pos[1] ** 2)
            return self._fin_attrs["fiducial_rmin"] < radius < self._fin_attrs["fiducial_rmax"]
        elif "fiducial_xmax" in self._fin_attrs:
            return (self._fin_attrs["fiducial_xmin"] < pos[0] < self._fin_attrs["fiducial_xmax"] and
                    self._fin_attrs["fiducial_ymin"] < pos[1] < self._fin_attrs["fiducial_ymax"])
        else:
            raise ValueError("Could not contruct fiducial volume from input file.")

    def _increase_signal(self, channel_id, factor):
        """
        increase the signal of a simulated station by a factor of x
        this is e.g. used to approximate a phased array concept with a single antenna

        Parameters
        ----------
        channel_id: int or None
            if None, all available channels will be modified
        """
        if channel_id is None:
            for electric_field in self._station.get_sim_station().get_electric_fields():
                electric_field.set_trace(electric_field.get_trace() * factor, sampling_rate=electric_field.get_sampling_rate())

        else:
            sim_channels = self._station.get_sim_station().get_electric_fields_for_channels([channel_id])
            for sim_channel in sim_channels:
                sim_channel.set_trace(sim_channel.get_trace() * factor, sampling_rate=sim_channel.get_sampling_rate())

    def _read_input_hdf5(self):
        """
        reads input file into memory
        """
        fin = h5py.File(self._inputfilename, 'r')
        self._fin = {}
        self._fin_stations = {}
        self._fin_attrs = {}
        for key, value in iteritems(fin):
            if isinstance(value, h5py._hl.group.Group):
                self._fin_stations[key] = {}
                for key2, value2 in iteritems(value):
                    self._fin_stations[key][key2] = np.array(value2)
            else:
                if len(value) and type(value[0]) == bytes:
                    self._fin[key] = np.array(value).astype('U')
                else:
                    self._fin[key] = np.array(value)
        for key, value in iteritems(fin.attrs):
            self._fin_attrs[key] = value

        fin.close()

    def _check_vertex_times(self):

        if 'vertex_times' in self._fin:
            return True
        else:
            warn_msg = 'The input file does not include vertex times. '
            warn_msg += 'Vertices from the same event will not be time-ordered.'
            logger.warning(warn_msg)
            return False

    def _calculate_signal_properties(self):
        if self._station.has_triggered():
            self._channelSignalReconstructor.run(self._evt, self._station, self._det)
            amplitudes = np.zeros(self._station.get_number_of_channels())
            amplitudes_envelope = np.zeros(self._station.get_number_of_channels())
            for channel in self._station.iter_channels():
                amplitudes[channel.get_id()] = channel.get_parameter(chp.maximum_amplitude)
                amplitudes_envelope[channel.get_id()] = channel.get_parameter(chp.maximum_amplitude_envelope)
            self._output_maximum_amplitudes[self._station.get_id()].append(amplitudes)
            self._output_maximum_amplitudes_envelope[self._station.get_id()].append(amplitudes_envelope)

    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=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=bool)
                sg['triggered'] = np.zeros(n_showers, dtype=bool)

    def _create_trigger_structures(self):

        if 'trigger_names' not in self._mout_attrs:
            self._mout_attrs['trigger_names'] = []
        extend_array = False
        for trigger in six.itervalues(self._station.get_triggers()):
            if trigger.get_name() not in self._mout_attrs['trigger_names']:
                self._mout_attrs['trigger_names'].append((trigger.get_name()))
                extend_array = True
        # the 'multiple_triggers' output array is not initialized in the constructor because the number of
        # 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=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=bool)
        elif extend_array:
            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
            # repeat for trigger times
            tmp_t = np.nan * np.zeros_like(tmp, dtype=float)
            tmp_t[:, 0:ny] = self._mout['trigger_times']
            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=bool)
#                 nx, ny = sg['multiple_triggers'].shape
#                 tmp[:, 0:ny] = sg['multiple_triggers']
#                 sg['multiple_triggers'] = tmp
        return extend_array

    def _save_triggers_to_hdf5(self, sg, local_shower_index, global_shower_index):

        extend_array = self._create_trigger_structures()
        # now we also need to create the trigger structure also in the sg (station group) dictionary that contains
        # 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=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=bool)
            nx, ny = sg['multiple_triggers'].shape
            tmp[:, 0:ny] = sg['multiple_triggers']
            sg['multiple_triggers'] = tmp
            # repeat for trigger times
            tmp_t = np.nan * np.zeros_like(tmp, dtype=float)
            tmp_t[:, :ny] = sg['trigger_times']
            sg['trigger_times'] = tmp_t

        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=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):
                multiple_triggers[iT] = self._station.get_trigger(trigger_name).has_triggered()
                trigger_times[iT] = self._station.get_trigger(trigger_name).get_trigger_time()
                for iSh in local_shower_index:  # now save trigger information per shower of the current station
                    sg['multiple_triggers'][iSh][iT] = self._station.get_trigger(trigger_name).has_triggered()
                    sg['trigger_times'][iSh][iT] = trigger_times[iT]
        for iSh, iSh2 in zip(local_shower_index, global_shower_index):  # now save trigger information per shower of the current station
            sg['triggered'][iSh] = np.any(sg['multiple_triggers'][iSh])
            self._mout['triggered'][iSh2] |= sg['triggered'][iSh]
            self._mout['multiple_triggers'][iSh2] |= sg['multiple_triggers'][iSh]
            self._mout['trigger_times'][iSh2] = np.fmin(
                self._mout['trigger_times'][iSh2], sg['trigger_times'][iSh])
        sg['event_id_per_shower'][local_shower_index] = self._evt.get_id()
        sg['event_group_id_per_shower'][local_shower_index] = self._evt.get_run_number()
        self._output_multiple_triggers_station[self._station_id].append(multiple_triggers)
        self._output_trigger_times_station[self._station_id].append(trigger_times)
        self._output_triggered_station[self._station_id].append(np.any(multiple_triggers))

    def get_Vrms(self):
        return self._Vrms

    def get_sampling_rate(self):
        return 1. / self._dt

    def get_bandwidth(self):
        return self._bandwidth

    def _check_if_was_pre_simulated(self):
        """
        checks if the same detector was simulated before (then we can save the ray tracing part)
        """
        self._was_pre_simulated = False

        if 'detector' in self._fin_attrs:
            if isinstance(self._det, detector.rnog_detector.Detector):
                if self._det.export_as_string() == self._fin_attrs['detector']:
                    self._was_pre_simulated = True
            else:
                with open(self._detectorfile, 'r') as fdet:
                    if fdet.read() == self._fin_attrs['detector']:
                        self._was_pre_simulated = True

        if self._was_pre_simulated:
            logger.status("the simulation was already performed with the same detector")

        return self._was_pre_simulated

    def _create_meta_output_datastructures(self):
        """
        creates the data structures of the parameters that will be saved into the hdf5 output file
        """
        self._mout = {}
        self._mout_attributes = {}
        self._mout['weights'] = np.zeros(self._n_showers)
        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 = {}
        self._output_triggered_station = {}
        self._output_event_group_ids = {}
        self._output_sub_event_ids = {}
        self._output_multiple_triggers_station = {}
        self._output_trigger_times_station = {}
        self._output_maximum_amplitudes = {}
        self._output_maximum_amplitudes_envelope = {}

        for station_id in self._station_ids:
            self._mout_groups[station_id] = {}
            sg = self._mout_groups[station_id]
            self._output_event_group_ids[station_id] = []
            self._output_sub_event_ids[station_id] = []
            self._output_triggered_station[station_id] = []
            self._output_multiple_triggers_station[station_id] = []
            self._output_trigger_times_station[station_id] = []
            self._output_maximum_amplitudes[station_id] = []
            self._output_maximum_amplitudes_envelope[station_id] = []

    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=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
        sg['event_group_id_per_shower'] = np.zeros(n_showers, dtype=int) * -1
        sg['launch_vectors'] = np.zeros((n_showers, n_antennas, nS, 3)) * np.nan
        sg['receive_vectors'] = np.zeros((n_showers, n_antennas, nS, 3)) * np.nan
        sg['polarization'] = np.zeros((n_showers, n_antennas, nS, 3)) * np.nan
        sg['travel_times'] = np.zeros((n_showers, n_antennas, nS)) * np.nan
        sg['travel_distances'] = np.zeros((n_showers, n_antennas, nS)) * np.nan
        if self._cfg['speedup']['amp_per_ray_solution']:
            sg['max_amp_shower_and_ray'] = np.zeros((n_showers, n_antennas, nS))
            sg['time_shower_and_ray'] = np.zeros((n_showers, n_antennas, nS))
        for parameter_entry in self._raytracer.get_output_parameters():
            if parameter_entry['ndim'] == 1:
                sg[parameter_entry['name']] = np.zeros((n_showers, n_antennas, nS)) * np.nan
            else:
                sg[parameter_entry['name']] = np.zeros((n_showers, n_antennas, nS, parameter_entry['ndim'])) * np.nan
        return sg

    def _read_input_particle_properties(self, idx=None):
        if idx is None:
            idx = self._primary_index
        self._event_group_id = self._fin['event_group_ids'][idx]

        self.input_particle = NuRadioReco.framework.particle.Particle(0)
        self.input_particle[simp.flavor] = self._fin['flavors'][idx]
        self.input_particle[simp.energy] = self._fin['energies'][idx]
        self.input_particle[simp.interaction_type] = self._fin['interaction_type'][idx]
        self.input_particle[simp.inelasticity] = self._fin['inelasticity'][idx]
        self.input_particle[simp.vertex] = np.array([self._fin['xx'][idx],
                                                     self._fin['yy'][idx],
                                                     self._fin['zz'][idx]])
        self.input_particle[simp.zenith] = self._fin['zeniths'][idx]
        self.input_particle[simp.azimuth] = self._fin['azimuths'][idx]
        self.input_particle[simp.inelasticity] = self._fin['inelasticity'][idx]
        self.input_particle[simp.n_interaction] = self._fin['n_interaction'][idx]
        if self._fin['n_interaction'][idx] <= 1:
            # parents before the neutrino and outgoing daughters without shower are currently not
            # simulated. The parent_id is therefore at the moment only rudimentarily populated.
            self.input_particle[simp.parent_id] = None  # primary does not have a parent

        self.input_particle[simp.vertex_time] = 0
        if 'vertex_times' in self._fin:
            self.input_particle[simp.vertex_time] = self._fin['vertex_times'][idx]

    def _read_input_shower_properties(self):
        """ read in the properties of the shower with index _shower_index from input """
        self._event_group_id = self._fin['event_group_ids'][self._shower_index]

        self._shower_vertex = np.array([self._fin['xx'][self._shower_index],
                                        self._fin['yy'][self._shower_index],
                                        self._fin['zz'][self._shower_index]])

        self._vertex_time = 0
        if 'vertex_times' in self._fin:
            self._vertex_time = self._fin['vertex_times'][self._shower_index]

    def _create_sim_station(self):
        """
        created an empyt sim_station object
        """
        # create NuRadioReco event structure
        self._sim_station = NuRadioReco.framework.sim_station.SimStation(self._station_id)
        self._sim_station.set_is_neutrino()

    def _create_sim_shower(self):
        """
        creates a sim_shower object and saves the meta arguments such as neutrino direction, shower energy and self.input_particle[flavor]
        """
        # create NuRadioReco event structure
        self._sim_shower = NuRadioReco.framework.radio_shower.RadioShower(self._shower_ids[self._shower_index])
        # save relevant neutrino properties
        self._sim_shower[shp.zenith] = self.input_particle[simp.zenith]
        self._sim_shower[shp.azimuth] = self.input_particle[simp.azimuth]
        self._sim_shower[shp.energy] = self._fin['shower_energies'][self._shower_index]
        self._sim_shower[shp.flavor] = self.input_particle[simp.flavor]
        self._sim_shower[shp.interaction_type] = self.input_particle[simp.interaction_type]
        self._sim_shower[shp.vertex] = self.input_particle[simp.vertex]
        self._sim_shower[shp.vertex_time] = self._vertex_time
        self._sim_shower[shp.type] = self._fin['shower_type'][self._shower_index]
        # TODO direct parent does not necessarily need to be the primary in general, but full
        # interaction chain is currently not populated in the input files.
        self._sim_shower[shp.parent_id] = self.primary.get_id()

    def _write_output_file(self, empty=False):
        folder = os.path.dirname(self._outputfilename)
        if not os.path.exists(folder) and folder != '':
            logger.warning(f"output folder {folder} does not exist, creating folder...")
            os.makedirs(folder)
        fout = h5py.File(self._outputfilename, 'w')

        if not empty:
            # here we add the first interaction to the saved events
            # if any of its children triggered

            # Careful! saved should be a copy of the triggered array, and not
            # a reference! saved indicates the interactions to be saved, while
            # triggered should indicate if an interaction has produced a trigger
            saved = np.copy(self._mout['triggered'])
            if 'n_interaction' in self._fin:  # if n_interactions is not specified, there are not parents
                parent_mask = self._fin['n_interaction'] == 1
                for event_id in np.unique(self._fin['event_group_ids']):
                    event_mask = self._fin['event_group_ids'] == event_id
                    if True in self._mout['triggered'][event_mask]:
                        saved[parent_mask & event_mask] = True

            logger.status("start saving events")
            # save data sets
            for (key, value) in iteritems(self._mout):
                fout[key] = value[saved]

            # save all data sets of the station groups
            for (key, value) in iteritems(self._mout_groups):
                sg = fout.create_group("station_{:d}".format(key))
                for (key2, value2) in iteritems(value):
                    sg[key2] = np.array(value2)[np.array(value['triggered'])]

            # save "per event" quantities
            if 'trigger_names' in self._mout_attrs:
                n_triggers = len(self._mout_attrs['trigger_names'])
                for station_id in self._mout_groups:
                    n_events_for_station = len(self._output_triggered_station[station_id])
                    if n_events_for_station > 0:
                        n_channels = self._det.get_number_of_channels(station_id)
                        sg = fout["station_{:d}".format(station_id)]
                        sg['event_group_ids'] = np.array(self._output_event_group_ids[station_id])
                        sg['event_ids'] = np.array(self._output_sub_event_ids[station_id])
                        sg['maximum_amplitudes'] = np.array(self._output_maximum_amplitudes[station_id])
                        sg['maximum_amplitudes_envelope'] = np.array(self._output_maximum_amplitudes_envelope[station_id])
                        sg['triggered_per_event'] = np.array(self._output_triggered_station[station_id])

                        # 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=bool)
                        for iE, values in enumerate(self._output_multiple_triggers_station[station_id]):
                            tmp[iE] = values
                        sg['multiple_triggers_per_event'] = tmp
                        tmp_t = np.nan * np.zeros_like(tmp, dtype=float)
                        for iE, values in enumerate(self._output_trigger_times_station[station_id]):
                            tmp_t[iE] = values
                        sg['trigger_times_per_event'] = tmp_t


        # save meta arguments
        for (key, value) in iteritems(self._mout_attrs):
            fout.attrs[key] = value

        if isinstance(self._det, detector.rnog_detector.Detector):
            fout.attrs['detector'] = self._det.export_as_string()
        else:
            with open(self._detectorfile, 'r') as fdet:
                fout.attrs['detector'] = fdet.read()

        if not empty:
            # save antenna position separately to hdf5 output
            for station_id in self._mout_groups:
                n_channels = self._det.get_number_of_channels(station_id)
                positions = np.zeros((n_channels, 3))
                for channel_id in range(n_channels):
                    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._integrated_channel_response[station_id].values())

            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)

        # save NuRadioMC and NuRadioReco versions
        from NuRadioReco.utilities import version
        import NuRadioMC
        fout.attrs['NuRadioMC_version'] = NuRadioMC.__version__
        fout.attrs['NuRadioMC_version_hash'] = version.get_NuRadioMC_commit_hash()

        if not empty:
            # now we also save all input parameters back into the out file
            for key in self._fin.keys():
                if key.startswith("station_"):
                    continue
                if not key in fout.keys():  # only save data sets that havn't been recomputed and saved already
                    if np.array(self._fin[key]).dtype.char == 'U':
                        fout[key] = np.array(self._fin[key], dtype=h5py.string_dtype(encoding='utf-8'))[saved]

                    else:
                        fout[key] = np.array(self._fin[key])[saved]

        for key in self._fin_attrs.keys():
            if not key in fout.attrs.keys():  # only save atrributes sets that havn't been recomputed and saved already
                if key not in ["trigger_names", "Tnoise", "Vrms", "bandwidth", "n_samples", "dt", "detector", "config"]:  # don't write trigger names from input to output file, this will lead to problems with incompatible trigger names when merging output files
                    fout.attrs[key] = self._fin_attrs[key]
        fout.close()

    def calculate_Veff(self):
        # calculate effective
        triggered = remove_duplicate_triggers(self._mout['triggered'], self._fin['event_group_ids'])
        n_triggered = np.sum(triggered)
        n_triggered_weighted = np.sum(self._mout['weights'][triggered])
        n_events = self._fin_attrs['n_events']
        logger.status(f'fraction of triggered events = {n_triggered:.0f}/{n_events:.0f} = {n_triggered / self._n_showers:.3f} (sum of weights = {n_triggered_weighted:.2f})')

        V = self._fin_attrs['volume']
        Veff = V * n_triggered_weighted / n_events
        logger.status(f"Veff = {Veff / units.km ** 3:.4g} km^3, Veffsr = {Veff * 4 * np.pi/units.km**3:.4g} km^3 sr")

    def _calculate_polarization_vector(self):
        """ calculates the polarization vector in spherical coordinates (eR, eTheta, ePhi)
        """
        if self._cfg['signal']['polarization'] == 'auto':
            polarization_direction = np.cross(self._launch_vector, np.cross(self._shower_axis, self._launch_vector))
            polarization_direction /= np.linalg.norm(polarization_direction)
            cs = cstrans.cstrafo(*hp.cartesian_to_spherical(*self._launch_vector))
            return cs.transform_from_ground_to_onsky(polarization_direction)
        elif self._cfg['signal']['polarization'] == 'custom':
            ePhi = float(self._cfg['signal']['ePhi'])
            eTheta = (1 - ePhi ** 2) ** 0.5
            v = np.array([0, eTheta, ePhi])
            return v / np.linalg.norm(v)
        else:
            msg = "{} for config.signal.polarization is not a valid option".format(self._cfg['signal']['polarization'])
            logger.error(msg)
            raise ValueError(msg)

    @property
    def _bandwidth_per_channel(self):
        warnings.warn("This variable `_bandwidth_per_channel` is deprecated. "
                      "Please use `integrated_channel_response` instead.", DeprecationWarning)
        return self._integrated_channel_response

    @_bandwidth_per_channel.setter
    def _bandwidth_per_channel(self, value):
        warnings.warn("This variable `_bandwidth_per_channel` is deprecated. "
                        "Please use `integrated_channel_response` instead.", DeprecationWarning)
        self._integrated_channel_response = value

    @property
    def integrated_channel_response(self):
        return self._integrated_channel_response

    @integrated_channel_response.setter
    def integrated_channel_response(self, value):
        self._integrated_channel_response = value
