Merge branch 'birefringence2' into refactor

NuRadioMC/test/SignalProp/T08test_emitter_birefringence.py

@@ -0,0 +1,74 @@
+import numpy as np
+from numpy import testing
+from NuRadioMC.SignalProp import analyticraytracing as ray
+from NuRadioMC.utilities import medium
+from NuRadioReco.utilities import units
+from NuRadioReco.framework import base_trace
+import NuRadioReco.framework.electric_field
+import logging
+from scipy.spatial.tests.test_qhull import points
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger('test_raytracing')
+
+import NuRadioReco.modules.io.eventReader as reader
+
+
+import matplotlib.pyplot as plt
+
+"""
+this unit test checks a full emitter simulation with birefringence
+"""
+
+#creating the reference traces
+"""
+event_reader = reader.eventReader()
+nur = 'emitter_sim_test/test_output.nur'
+event_reader.begin([nur])
+list_traces = []
+
+for event in event_reader.run():
+
+    station = event.get_station(51)
+    list_trace = []
+
+    for channel in station.iter_channels():
+
+        trace = channel.get_trace()
+        list_trace.append(trace)
+
+    event_traces = np.vstack(list_trace)
+    list_traces.append(event_traces)
+
+traces = np.dstack(list_traces)
+traces = np.transpose(traces, (2, 0, 1))
+
+np.save('reference_emitter.npy', traces)
+"""
+
+event_reader = reader.eventReader()
+nur = 'emitter_sim_test/test_output.nur'
+event_reader.begin([nur])
+list_traces = []
+
+for event in event_reader.run():
+
+    station = event.get_station(51)
+    list_trace = []
+
+    for channel in station.iter_channels():
+
+        trace = channel.get_trace()
+        list_trace.append(trace)
+
+    event_traces = np.vstack(list_trace)
+    list_traces.append(event_traces)
+
+traces = np.dstack(list_traces)
+traces = np.transpose(traces, (2, 0, 1))
+
+reference_array = np.load('reference_emitter.npy')
+
+testing.assert_allclose(traces, reference_array, atol=5e-5  * units.V / units.m, rtol=1e-7)
+print('T08test_emitter_birefringence passed without issues')
+
+

NuRadioMC/test/SignalProp/emitter_sim_test/T01_sim_events.py

@@ -0,0 +1,56 @@
+from __future__ import absolute_import, division, print_function
+import numpy as np
+from NuRadioReco.utilities import units
+from six import iterkeys, iteritems
+from scipy import constants
+from scipy.integrate import quad
+from scipy.interpolate import interp1d
+from scipy.optimize import fsolve
+import h5py
+from NuRadioMC.EvtGen.generator import write_events_to_hdf5
+import logging
+logger = logging.getLogger("EventGen")
+logging.basicConfig()
+
+VERSION_MAJOR = 1
+VERSION_MINOR = 1
+
+
+def generate_my_events(filename, n_events):
+    """
+    Event generator skeleton
+
+    Parameters
+    ----------
+    filename: string
+        the output filename of the hdf5 file
+    n_events: int
+        number of events to generate
+    """
+
+    # first set the meta attributes
+    attributes = {}
+    n_events = int(n_events)
+    attributes['n_events'] = n_events
+    attributes['start_event_id'] = 0
+    attributes['simulation_mode'] = 'emitter'
+
+    # now generate the events and fill all required data sets
+    # here we fill all data sets with dummy values
+    data_sets = {}
+
+    # define the emitter positions. X/Y are the easting/northing coordinates of the SPICE core (https://journals.aps.org/prd/pdf/10.1103/PhysRevD.105.123012)
+    data_sets["xx"] = np.ones(n_events) * 12911 * units.m
+    data_sets["yy"] = np.ones(n_events) * 14927.3 * units.m
+    data_sets["zz"] = -np.linspace(0, 2000, n_events) * units.m
+
+    data_sets["event_group_ids"] = np.arange(n_events)
+    data_sets["shower_ids"] = np.arange(n_events)
+
+    data_sets['emitter_model'] = np.array(['efield_idl1_spice'] * n_events)
+    data_sets['emitter_amplitudes'] = np.ones(n_events) # for the efield_idl1_spice model, this quantity is the relative rescaling of the measured amplitudes
+
+    write_events_to_hdf5(filename, data_sets, attributes)
+
+if __name__ == "__main__":
+    generate_my_events("test_input.hdf5", 10)

NuRadioMC/test/SignalProp/emitter_sim_test/T02_run_simulation.py

@@ -0,0 +1,53 @@
+from __future__ import absolute_import, division, print_function
+import argparse
+import datetime
+# import detector simulation modules
+import NuRadioReco.modules.efieldToVoltageConverter
+import NuRadioReco.modules.channelResampler
+import NuRadioReco.modules.channelGenericNoiseAdder
+import NuRadioReco.modules.triggerTimeAdjuster
+import NuRadioReco.modules.ARIANNA.hardwareResponseIncorporator
+import NuRadioReco.modules.trigger.highLowThreshold
+from NuRadioReco.utilities import units
+from NuRadioMC.simulation import simulation
+import logging
+logging.basicConfig(level=logging.WARNING)
+logger = logging.getLogger("runstrawman")
+
+# initialize detector sim modules
+efieldToVoltageConverter = NuRadioReco.modules.efieldToVoltageConverter.efieldToVoltageConverter()
+efieldToVoltageConverter.begin()
+channelResampler = NuRadioReco.modules.channelResampler.channelResampler()
+channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder()
+hardwareResponseIncorporator = NuRadioReco.modules.ARIANNA.hardwareResponseIncorporator.hardwareResponseIncorporator()
+triggerSimulator = NuRadioReco.modules.trigger.highLowThreshold.triggerSimulator()
+triggerTimeAdjuster = NuRadioReco.modules.triggerTimeAdjuster.triggerTimeAdjuster()
+triggerSimulator.begin(log_level=logging.WARNING)
+
+class mySimulation(simulation.simulation):
+
+    def _detector_simulation_filter_amp(self, evt, station, det):
+        hardwareResponseIncorporator.run(evt, station, det, sim_to_data=True)
+
+    def _detector_simulation_trigger(self, evt, station, det):
+        triggerSimulator.run(evt, station, det,
+                           threshold_high=0.5 * units.mV,
+                           threshold_low=-0.5 * units.mV,
+                           high_low_window=5 * units.ns,
+                           coinc_window=30 * units.ns,
+                           number_concidences=2,
+                           triggered_channels=range(8))
+        triggerTimeAdjuster.run(evt, station, det)
+
+if __name__ == "__main__":
+    sim = mySimulation(inputfilename='test_input.hdf5',
+                                outputfilename='test_output.hdf5',
+                                detectorfile='test_detector.json',
+                                outputfilenameNuRadioReco='test_output.nur',
+                                config_file='test_config.yaml',
+                                #log_level=logging.WARNING,
+                                log_level=logging.ERROR,
+                                evt_time=datetime.datetime(2018, 12, 30),
+                                file_overwrite=True)
+    sim.run()
+

NuRadioMC/test/SignalProp/emitter_sim_test/test_config.yaml

@@ -0,0 +1,25 @@
+noise: False  # specify if simulation should be run with or without noise
+sampling_rate: 5.  # sampling rate in GHz used internally in the simulation. At the end the waveforms will be downsampled to the sampling rate specified in the detector description
+
+speedup:
+  minimum_weight_cut: 0 # disable neutrino weight cut -> we're not simulating neutrinos here
+  delta_C_cut: 10  # more than 360deg, disable cone angle cut
+  min_efield_amplitude: 0
+  distance_cut: False
+
+propagation:
+  module: analytic
+  ice_model: southpole_2015
+  attenuation_model: SP1
+  attenuate_ice: True # if True apply the frequency dependent attenuation due to propagating through ice. (Note: The 1/R amplitude scaling will be applied in either case.)
+  n_freq: 25  # the number of frequencies where the attenuation length is calculated for. The remaining frequencies will be determined from a linear interpolation between the reference frequencies. The reference frequencies are equally spaced over the complet frequency range. 
+  birefringence: True
+  birefringence_propagation: 'analytical'
+  birefringence_model: 'southpole_A'
+  angle_to_iceflow: -131
+
+signal:
+  model: spherical
+  zerosignal: False  # if True, the signal is set to zero. This is useful to study 'noise' only simulations
+  polarization: custom # can be either 'auto' or 'custom'
+  ePhi: 0.  # only used if 'polarization = custom', fraction of ePhi component, the eTheta component is eTheta = (1 - ePhi**2)**0.5