diff --git a/.github/workflows/run_tests.yaml b/.github/workflows/run_tests.yaml index 9d092e45c..28ce3711e 100644 --- a/.github/workflows/run_tests.yaml +++ b/.github/workflows/run_tests.yaml @@ -195,4 +195,5 @@ jobs: run: | export GSLDIR=$(gsl-config --prefix) export PYTHONPATH=$(pwd):$PYTHONPATH + python install_dev.py --dev cr_interpolator --no-interactive NuRadioReco/test/test_examples.sh diff --git a/NuRadioReco/examples/AirShowerFluenceXmaxReco/A01createDetectorJson.py b/NuRadioReco/examples/AirShowerFluenceXmaxReco/A01createDetectorJson.py new file mode 100644 index 000000000..a95ca592f --- /dev/null +++ b/NuRadioReco/examples/AirShowerFluenceXmaxReco/A01createDetectorJson.py @@ -0,0 +1,101 @@ +""" +Helper script that generates the detector description for +dual-polarized antennas on a square grid. This is a generic +example for typical air-shower arrays such as AERA, Auger-Radio, +LOFAR, TUNKA-REX, GRAND, etc. +""" + +import json +import copy +import numpy as np +from radiotools import helper as hp +from NuRadioReco.utilities import units + +class NpEncoder(json.JSONEncoder): + + def default(self, obj): + if isinstance(obj, np.integer): + return int(obj) + elif isinstance(obj, np.floating): + return float(obj) + elif isinstance(obj, np.ndarray): + return obj.tolist() + else: + return super(NpEncoder, self).default(obj) + + +antenna = "SKALA_InfFirn" +ori_EW = [0,0,90*units.deg, 90*units.deg] +ori_NW = [0,0,90*units.deg, 0] + +# antenna = "LOFAR_LBA" +# ori_EW = [90 * units.deg, 135 * units.deg, 0, 0] +# ori_NW = [90 * units.deg, (135+180) * units.deg, 0, 0] + +orientation_theta, orientation_phi, rotation_theta, rotation_phi = ori_NW +a1 = hp.spherical_to_cartesian(orientation_theta, orientation_phi) +a2 = hp.spherical_to_cartesian(rotation_theta, rotation_phi) +a3 = np.cross(a1, a2) +if np.linalg.norm(a3) < 0.9: + raise AssertionError("the two vectors that define the antenna orientation are not othorgonal to each other") + +ctmp = { + "adc_n_samples": 512, + "adc_sampling_frequency": 1, + "channel_id": 0, + "commission_time": "{TinyDate}:2017-11-01T00:00:00", + "decommission_time": "{TinyDate}:2038-01-01T00:00:00", +} + +stmp = { + "commission_time": "{TinyDate}:2017-11-04T00:00:00", + "decommission_time": "{TinyDate}:2038-01-01T00:00:00", + "pos_altitude": 0, + "pos_site": "LOFAR", +} + +det = {} +det['channels'] = {} +det['stations'] = {} + +xx = np.arange(-0.15,0.151,0.05) * units.km +yy = np.arange(-0.15,0.151,0.05) * units.km +# yy = [0] +z = 1 * units.cm + +station_id = 1 +channel_id_counter = 0 +channel_group_counter = 0 +for x in xx: + for y in yy: + channel_group_counter += 1 + for ori in [ori_EW, ori_NW]: + channel_id_counter += 1 + ori_theta, ori_phi, rot_theta, rot_phi = ori + cab_length = 0 + cab_time_delay = 0 + Tnoise = 300 + channel = copy.copy(ctmp) + channel['ant_position_x'] = x + channel['ant_position_y'] = y + channel['ant_position_z'] = z + channel['ant_type'] = antenna + channel['channel_id'] = channel_id_counter + channel['channel_group_id'] = channel_group_counter + channel["ant_orientation_phi"] = ori_phi/units.deg + channel["ant_orientation_theta"] = ori_theta/units.deg + channel["ant_rotation_phi"] = rot_phi/units.deg + channel["ant_rotation_theta"] = rot_theta/units.deg + channel["cab_length"] = cab_length + channel["cab_time_delay"] = cab_time_delay + channel["noise_temperature"] = Tnoise + channel["station_id"] = station_id + det['channels']["{:d}".format(channel_id_counter)] = channel + det['stations'][f"{station_id:d}"] = copy.copy(stmp) + det['stations'][f"{station_id:d}"]['station_id'] = station_id + det['stations'][f"{station_id:d}"]['pos_easting'] = 0 + det['stations'][f"{station_id:d}"]['pos_northing'] = 0 + + with open(f"grid_array_{antenna}.json", 'w') as fout: + json.dump(det, fout, indent=4, sort_keys=True, cls=NpEncoder) + diff --git a/NuRadioReco/examples/AirShowerFluenceXmaxReco/A02CoREASFluenceXmaxReco.py b/NuRadioReco/examples/AirShowerFluenceXmaxReco/A02CoREASFluenceXmaxReco.py new file mode 100644 index 000000000..3cade9076 --- /dev/null +++ b/NuRadioReco/examples/AirShowerFluenceXmaxReco/A02CoREASFluenceXmaxReco.py @@ -0,0 +1,289 @@ +""" +Reconstruction of the position of the position of the shower maximum Xmax +using the "LOFAR"-style analysis where the best CoREAS simulation is found +which matches the data best by comparing the energy fluence footprint. + +The procedure is the following (following the description in S. Buitink Phys. Rev. D 90, 082003 https://doi.org/10.1103/PhysRevD.90.082003) + * One star-shape CoREAS simulation is read in and interpolate to obtain the electric field + at every antenna position. This is handled by the `readCoREASDetector` module + * A full detector simulation is performed, i.e., simulation of antenna and signal chain + response, adding of noise, electric field reconstruction (using standard unfolding), and + calculation of the energy fluence at every antenna position + * we loop over a set of CoREAS simulations (excluding the simulation we used to generate the data) + of the same air-shower direction and energy and + fit the fluence footprint to the simulated data to determine the core position and amplitude + normalization factor. The chi^2 of the best fit is saved. + * the best fit chi^2 values are plotted against the true Xmax value from the simulations + * a parabola is fitted to determine the most likely Xmax value +""" + +import os +import glob +import logging +import numpy as np +from matplotlib import pyplot as plt +from matplotlib import cm +from scipy import constants +from scipy import optimize as opt +from astropy import time +from NuRadioReco.detector import detector +from NuRadioReco.modules.io.coreas import coreas +from NuRadioReco.utilities import units +from NuRadioReco.utilities.dataservers import download_from_dataserver +import NuRadioReco.modules.io.coreas.readCoREASDetector +import NuRadioReco.modules.io.coreas.simulationSelector +import NuRadioReco.modules.efieldToVoltageConverter +import NuRadioReco.modules.ARIANNA.hardwareResponseIncorporator +import NuRadioReco.modules.channelGenericNoiseAdder +import NuRadioReco.modules.trigger.simpleThreshold +import NuRadioReco.modules.channelBandPassFilter +import NuRadioReco.modules.eventTypeIdentifier +import NuRadioReco.modules.channelStopFilter +import NuRadioReco.modules.channelSignalReconstructor +import NuRadioReco.modules.correlationDirectionFitter +import NuRadioReco.modules.voltageToEfieldConverter +import NuRadioReco.modules.electricFieldSignalReconstructor +import NuRadioReco.modules.electricFieldBandPassFilter +import NuRadioReco.modules.voltageToAnalyticEfieldConverter +import NuRadioReco.modules.voltageToEfieldConverterPerChannelGroup +import NuRadioReco.modules.channelResampler +import NuRadioReco.modules.electricFieldResampler +import NuRadioReco.modules.io.eventWriter +from NuRadioReco.detector import antennapattern +from NuRadioReco.framework.parameters import electricFieldParameters as efp +from NuRadioReco.framework.parameters import showerParameters as shp + + +if not os.path.exists("plots"): + os.makedirs("plots") + +# download coreas simulations +remote_path = "data/CoREAS/LOFAR/evt_00001" +path = 'data/CoREAS/LOFAR/evt_00001' +for i in [0,1,2,3,4,5,6,7,9,10,11,12,13,14,15,16,17,18,19,21,22,23,24,25,26,27,28,29]: + filename = os.path.join(remote_path, f"SIM0000{i:02d}.hdf5") + if not os.path.exists(filename): + download_from_dataserver(os.path.join(remote_path, f"SIM0000{i:02d}.hdf5"), + os.path.join(path, f"SIM0000{i:02d}.hdf5")) + +# Initialize detector +det = detector.Detector(json_filename="grid_array_SKALA_InfFirn.json", + assume_inf=False, antenna_by_depth=False) +det.update(time.Time("2023-01-01T00:00:00", format='isot', scale='utc')) + + +# Initialize all modules that are used in the reconstruction pipeline +electricFieldSignalReconstructor = NuRadioReco.modules.electricFieldSignalReconstructor.electricFieldSignalReconstructor() +electricFieldSignalReconstructor.begin() +simulationSelector = NuRadioReco.modules.io.coreas.simulationSelector.simulationSelector() +simulationSelector.begin() +efieldToVoltageConverter = NuRadioReco.modules.efieldToVoltageConverter.efieldToVoltageConverter(log_level=logging.INFO) +efieldToVoltageConverter.begin(debug=False, pre_pulse_time=0, post_pulse_time=0) +hardwareResponseIncorporator = NuRadioReco.modules.ARIANNA.hardwareResponseIncorporator.hardwareResponseIncorporator() +channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() +triggerSimulator = NuRadioReco.modules.trigger.simpleThreshold.triggerSimulator() +triggerSimulator.begin() +channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() +channelBandPassFilter.begin() +eventTypeIdentifier = NuRadioReco.modules.eventTypeIdentifier.eventTypeIdentifier() +channelSignalReconstructor = NuRadioReco.modules.channelSignalReconstructor.channelSignalReconstructor() +channelSignalReconstructor.begin() +voltageToEfieldConverter = NuRadioReco.modules.voltageToEfieldConverter.voltageToEfieldConverter() +electricFieldBandPassFilter = NuRadioReco.modules.electricFieldBandPassFilter.electricFieldBandPassFilter() +voltageToEfieldConverterPerChannelGroup = NuRadioReco.modules.voltageToEfieldConverterPerChannelGroup.voltageToEfieldConverterPerChannelGroup() +voltageToEfieldConverterPerChannelGroup.begin(use_MC_direction=True) + +coreas_reader = NuRadioReco.modules.io.coreas.readCoREASDetector.readCoREASDetector() + +antenna_pattern_provider = antennapattern.AntennaPatternProvider() + +core = np.array([0, 0, 0]) * units.m + + +filt_settings = {'passband': [80 * units.MHz, 300 * units.MHz], + 'filter_type': 'butter', + 'order': 10} + +input_file = os.path.join(path, "SIM000000.hdf5") +coreas_reader.begin(input_file) +# we only need a single realization of the shower, so we set the core position to zero for simplicity +for iE, evt in enumerate(coreas_reader.run(det, [core])): + station = evt.get_station() + simsh = evt.get_first_sim_shower() + + if 0: # example plotting of the interpolated efields. + ss = station.get_sim_station() + for efield in ss.get_electric_fields(): + fig, ax = plt.subplots(1, 1) + ax.plot(efield.get_times(), efield.get_trace()[1], label='eTheta') + ax.plot(efield.get_times(), efield.get_trace()[2], label='ePhi') + ax.set_xlabel('time [ns]') + ax.set_ylabel('electric field [V/m]') + ax.legend() + pos = det.get_relative_position(station.get_id(), efield.get_unique_identifier()[0][0]) + ax.set_title(f"E = {simsh[shp.energy]:.2g}eV, obs = {simsh[shp.observation_level]:.0f}m, pos = ({pos[0]:.0f}m , {pos[1]:.0f}m, {pos[2]:.0f}m) zenith = {efield[efp.zenith]/units.deg:.0f}deg, azimuth = {efield[efp.azimuth]/units.deg:.0f}deg", fontsize="small") + fig.tight_layout() + # fig.savefig(f"plots/efield_{efield.get_unique_identifier()[0][0]}.png") + plt.show() + + # apply antenna response + efieldToVoltageConverter.run(evt, station, det) + + # approximate the rest of the signal chain with a bandpass filter + channelBandPassFilter.run(evt, station, det, **filt_settings) + + # add thermal noise of fixed noise temperature + Tnoise = 300 # in Kelvin + + # calculate Vrms and normalize such that after filtering the correct Vrms is obtained + min_freq = 0 + max_freq = 0.5 * det.get_sampling_frequency(station.get_id(), 1) + ff = np.linspace(0, max_freq, 10000) + filt = channelBandPassFilter.get_filter(ff, station.get_id(), None, det, **filt_settings) + bandwidth = np.trapz(np.abs(filt) ** 2, ff) + Vrms = (Tnoise * 50 * constants.k * bandwidth / units.Hz) ** 0.5 + amplitude = Vrms/ (bandwidth / max_freq) ** 0.5 + channelGenericNoiseAdder.run(evt, station, det, type="rayleigh", amplitude=1 * amplitude, + min_freq=min_freq, max_freq=max_freq) + + # simulate if the air shower triggers the station. Requireing a 10 sigma simple threshold trigger + triggerSimulator.run(evt, station, det, number_concidences=1, threshold=10*Vrms) + if station.get_trigger('default_simple_threshold').has_triggered(): + + eventTypeIdentifier.run(evt, station, "forced", 'cosmic_ray') + channelSignalReconstructor.run(evt, station, det) + + # reconstruct the electric field for each dual-polarized antenna through standard unfolding + voltageToEfieldConverterPerChannelGroup.run(evt, station, det) + simsh = evt.get_sim_shower(0) + + # calcualte the energy fluence from every reconstructed efield and save in arrays + electricFieldSignalReconstructor.run(evt, station, det) + ff = [] + ff_error = [] + pos = [] + for efield in station.get_electric_fields(): + ff.append(np.sum(efield[efp.signal_energy_fluence])) + ff_error.append(np.sum(efield.get_parameter_error(efp.signal_energy_fluence))) + pos.append(efield.get_position()) + pos = np.array(pos) + ff = np.array(ff) + ff_error = np.array(ff_error) + ff_error[ff_error == 0] = ff[ff_error != 0].min() # ensure that the error is always larger than 0 to avoid division by zero in objective function + + if 1: # plot energy fluence footprint + fig, ax = plt.subplots(1, 1) + sc = ax.scatter(pos[:, 0], pos[:, 1], c=ff, cmap=cm.gnuplot2_r, marker='o', edgecolors='k') + cbar = fig.colorbar(sc, ax=ax, orientation='vertical') + cbar.set_label('fluence [eV/m^2]') + ax.set_title(f"E = {simsh[shp.energy]:.2g}eV, obs = {simsh[shp.observation_level]:.0f}m, core = {core[0]:.0f}m, {core[1]:.0f}m, {core[2]:.0f}m") + ax.set_xlabel('x [m]') + ax.set_ylabel('y [m]') + fig.savefig(f"plots/footprint_{iE}_core{core[0]:0f}_{core[1]:0f}.png") + plt.close("all") + # plt.show() + + if 0: # plot individual reconstructed efields + for efield in station.get_electric_fields(): + if np.abs(efield.get_position()[1]) < 1: + fig, ax = plt.subplots(1, 1) + ax.plot(efield.get_times(), efield.get_trace()[1], label='eTheta') + ax.plot(efield.get_times(), efield.get_trace()[2], label='ePhi') + ax.set_xlabel('time [ns]') + ax.set_ylabel('electric field [V/m]') + ax.legend() + tpos = det.get_relative_position(station.get_id(), efield.get_unique_identifier()[0][0]) + ax.set_title(f"E = {simsh[shp.energy]:.2g}eV, obs = {simsh[shp.observation_level]:.0f}m, pos = ({tpos[0]:.0f}m , {tpos[1]:.0f}m, {tpos[2]:.0f}m) zenith = {efield[efp.zenith]/units.deg:.0f}deg, azimuth = {efield[efp.azimuth]/units.deg:.0f}deg", fontsize="small") + fig.tight_layout() + fig.savefig(f"plots/efield/efield_core{iE}_{efield.get_unique_identifier()[0][0]}.png") + plt.close("all") + plt.show() + + # perform Xmax fit: + + input_files = glob.glob(os.path.join(path, "*.hdf5")) + N_showers = len(input_files) + + # initialize arrays to store the fit results + fmin = np.zeros(N_showers) + success = np.zeros(N_showers, dtype=bool) + Xmax = np.zeros(N_showers) + # loop through all + for i, filename2 in enumerate(input_files): + # skip the CoREAS simulation that was used to generate the event under test + if(filename2 == input_file): + continue + + # read in CoREAS simulation + evt2 = coreas.read_CORSIKA7(filename2) + Xmax[i] = evt2.get_first_sim_shower()[shp.shower_maximum] + # preprocess simulation, apply the same bandpass filter + # that the fake data has. This is required because the energy + # fluence depends on the bandwidth/choosen filter + ss = evt2.get_station(0).get_sim_station() + electricFieldBandPassFilter.run(evt2, ss, det, **filt_settings) + electricFieldSignalReconstructor.run(evt2, ss, None) + # Initialize interpolator for core position fit + coreasInterpolator = coreas.coreasInterpolator(evt2) + coreasInterpolator.initialize_fluence_interpolator() + + # define objective function for core position and amplitude fit + def obj(xyA): + A = xyA[2] # normalization + core = np.array([xyA[0], xyA[1], 0]) * units.m # core position + ff_true = np.zeros_like(ff) + for ip, p in enumerate(pos): + ff_true[ip] = coreasInterpolator.get_interp_fluence_value(p, core) + ff_true *= A + return np.sum((ff - ff_true) ** 2/ff_error**2) + + # perform minimization + res = opt.minimize(obj, [core[0], core[1], 1], method='Nelder-Mead') + # save fit results + fmin[i] = res.fun + success[i] = res.success + print(f"fitting footprint {i} from file {filename2}") + print(res) + + + # End loop CoREAS simulations + # now we can plot the best chi2 vs. Xmax + ndf = len(ff) - 3 + m = fmin > 0 # let's exclude unsuccessful footprint fits + xx = Xmax[m]/units.g * units.cm**2 + yy = fmin[m]/ndf + + # define parabola in standard-normal form + def func(x, a, b, c): + return a * (x-b)**2 + np.abs(c) + + # define objective function + def obj(p, xx, yy): + return np.sum((yy - func(xx, *p))**2) + + # fit parabola + xmax_best_fit = xx[yy.argmin()] + mask = (xx > (xmax_best_fit - 100)) & (xx < (xmax_best_fit + 100)) & success[m] + res = opt.minimize(obj, [1, xmax_best_fit, 0], args=(xx[mask], yy[mask]), method='Nelder-Mead') + popt = res.x + + xmax_best_fit = popt[1] + true_xmax = simsh[shp.shower_maximum]/units.g * units.cm**2 + xmax_diff = np.abs(xmax_best_fit - true_xmax) + + print(f"true Xmax = {true_xmax:.2f} g/cm^2, best fit Xmax = {xmax_best_fit:.2f} g/cm^2, diff = {xmax_diff:.2f} g/cm^2") + + if 1: # plot parabola + fig, ax = plt.subplots(1, 1) + ax.plot(xx, yy, 'o') + ax.plot(xx[~success[m]], yy[~success[m]], 'ro') + xxx = np.linspace(xx[mask].min(), xx[mask].max(), 1000) + ax.plot(xxx, func(xxx, *popt)) + ax.vlines(simsh[shp.shower_maximum]/units.g * units.cm**2, 0, yy.max(), linestyles='dashed', label='true Xmax') + ax.set_xlabel('Xmax [g/cm^2]') + ax.set_ylabel('chi2 / ndf') + ax.set_title(f"true Xmax = {true_xmax:.2f} g/cm^2, best fit Xmax = {xmax_best_fit:.2f} g/cm^2, diff = {xmax_diff:.2f} g/cm^2", fontsize="small") + fig.tight_layout() + fig.savefig(f"plots/Xmax_fit_{iE}.png") + plt.close("all") + # plt.show() \ No newline at end of file diff --git a/NuRadioReco/examples/AirShowerFluenceXmaxReco/README.md b/NuRadioReco/examples/AirShowerFluenceXmaxReco/README.md new file mode 100644 index 000000000..88640bc52 --- /dev/null +++ b/NuRadioReco/examples/AirShowerFluenceXmaxReco/README.md @@ -0,0 +1,20 @@ +Reconstruction of the position of the position of the shower maximum Xmax +using the "LOFAR"-style analysis where the best CoREAS simulation is found +which matches the data best by comparing the energy fluence footprint. + +The procedure is the following (following the description in S. Buitink Phys. Rev. D 90, 082003 https://doi.org/10.1103/PhysRevD.90.082003) + * One star-shape CoREAS simulation is read in and interpolate to obtain the electric field + at every antenna position. This is handled by the `readCoREASDetector` module + * A full detector simulation is performed, i.e., simulation of antenna and signal chain + response, adding of noise, electric field reconstruction (using standard unfolding), and + calculation of the energy fluence at every antenna position + * we loop over a set of CoREAS simulations (excluding the simulation we used to generate the data) + of the same air-shower direction and energy and + fit the fluence footprint to the simulated data to determine the core position and amplitude + normalization factor. The chi^2 of the best fit is saved. + * the best fit chi^2 values are plotted against the true Xmax value from the simulations + * a parabola is fitted to determine the most likely Xmax value + + + A01createDetectorJson.py creastes the detector description + A02CoREASFluenceXmaxReco.py downloads example CORSIKA7 simulations (verticla shower at the LOFAR site) and performs the fit \ No newline at end of file diff --git a/NuRadioReco/examples/FullReconstruction.py b/NuRadioReco/examples/FullReconstruction.py deleted file mode 100644 index 73334d16b..000000000 --- a/NuRadioReco/examples/FullReconstruction.py +++ /dev/null @@ -1,162 +0,0 @@ -import os -import sys -import datetime -import matplotlib.pyplot as plt -import astropy -from NuRadioReco.utilities import units -from NuRadioReco.detector import detector -import NuRadioReco.modules.io.coreas.readCoREAS -import NuRadioReco.modules.io.coreas.simulationSelector -import NuRadioReco.modules.efieldToVoltageConverter -import NuRadioReco.modules.ARIANNA.hardwareResponseIncorporator -import NuRadioReco.modules.channelGenericNoiseAdder -import NuRadioReco.modules.trigger.simpleThreshold -import NuRadioReco.modules.channelBandPassFilter -import NuRadioReco.modules.eventTypeIdentifier -import NuRadioReco.modules.channelStopFilter -import NuRadioReco.modules.channelSignalReconstructor -import NuRadioReco.modules.correlationDirectionFitter -import NuRadioReco.modules.voltageToEfieldConverter -import NuRadioReco.modules.electricFieldSignalReconstructor -import NuRadioReco.modules.voltageToAnalyticEfieldConverter -import NuRadioReco.modules.channelResampler -import NuRadioReco.modules.electricFieldResampler -import NuRadioReco.modules.io.eventWriter - -# Logging -import logging -logger = logging.getLogger("NuRadioReco.FullReconstruction") # Logging level is globally controlled - -plt.switch_backend('agg') - - -""" -Here, we show an example reconstruction of CoREAS data. A variety of modules -are being used. Please refer to details in the modules themselves. - -Input parameters (all with a default provided) ---------------------- - -Command line input: - python FullReconstruction.py station_id input_file detector_file - -station_id: int - station id to be used, default 32 -input_file: str - CoREAS simulation file, default example data -detector_file: str - path to json detector database, default given -""" - -try: - station_id = int(sys.argv[1]) # specify station id - input_file = sys.argv[2] # file with coreas simulations -except: - print("Usage: python FullReconstruction.py station_id input_file detector") - station_id = 32 - input_file = "example_data/example_event.h5" - print("Using default station {}".format(32)) - -if(station_id == 32): - triggered_channels = [0, 1, 2, 3] - used_channels_efield = [0, 1, 2, 3] - used_channels_fit = [0, 1, 2, 3] - channel_pairs = ((0, 2), (1, 3)) -else: - print("Default channels not defined for station_id != 32") - -try: - detector_file = sys.argv[3] - print("Using {0} as detector".format(detector_file)) -except: - print("Using default file for detector") - detector_file = '../examples/example_data/arianna_station_32.json' - - -det = detector.Detector(json_filename=detector_file) # detector file -det.update(datetime.datetime(2018, 10, 1)) - -dir_path = os.path.dirname(os.path.realpath(__file__)) # get the directory of this file - -# initialize all modules that are needed for processing -# provide input parameters that are to remain constant during processung -readCoREAS = NuRadioReco.modules.io.coreas.readCoREAS.readCoREAS() -readCoREAS.begin([input_file], station_id, n_cores=10, max_distance=None) -simulationSelector = NuRadioReco.modules.io.coreas.simulationSelector.simulationSelector() -simulationSelector.begin() -efieldToVoltageConverter = NuRadioReco.modules.efieldToVoltageConverter.efieldToVoltageConverter() -efieldToVoltageConverter.begin(debug=False) -hardwareResponseIncorporator = NuRadioReco.modules.ARIANNA.hardwareResponseIncorporator.hardwareResponseIncorporator() -channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() -triggerSimulator = NuRadioReco.modules.trigger.simpleThreshold.triggerSimulator() -triggerSimulator.begin() -channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() -channelBandPassFilter.begin() -eventTypeIdentifier = NuRadioReco.modules.eventTypeIdentifier.eventTypeIdentifier() -channelStopFilter = NuRadioReco.modules.channelStopFilter.channelStopFilter() -channelSignalReconstructor = NuRadioReco.modules.channelSignalReconstructor.channelSignalReconstructor() -channelSignalReconstructor.begin() -correlationDirectionFitter = NuRadioReco.modules.correlationDirectionFitter.correlationDirectionFitter() -voltageToEfieldConverter = NuRadioReco.modules.voltageToEfieldConverter.voltageToEfieldConverter() - -electricFieldSignalReconstructor = NuRadioReco.modules.electricFieldSignalReconstructor.electricFieldSignalReconstructor() -electricFieldSignalReconstructor.begin() - -voltageToAnalyticEfieldConverter = NuRadioReco.modules.voltageToAnalyticEfieldConverter.voltageToAnalyticEfieldConverter() -voltageToAnalyticEfieldConverter.begin() - -electricFieldResampler = NuRadioReco.modules.electricFieldResampler.electricFieldResampler() -electricFieldResampler.begin() - -channelResampler = NuRadioReco.modules.channelResampler.channelResampler() -channelResampler.begin() -eventWriter = NuRadioReco.modules.io.eventWriter.eventWriter() -output_filename = "MC_example_station_{}.nur".format(station_id) -eventWriter.begin(output_filename) - -# Loop over all events in file as initialized in readCoRREAS and perform analysis -for iE, evt in enumerate(readCoREAS.run(detector=det)): - logger.info("processing event {:d} with id {:d}".format(iE, evt.get_id())) - station = evt.get_station(station_id) - station.set_station_time(astropy.time.Time('2019-01-01T00:00:00')) - det.update(station.get_station_time()) - if simulationSelector.run(evt, station.get_sim_station(), det): - - efieldToVoltageConverter.run(evt, station, det) - - hardwareResponseIncorporator.run(evt, station, det, sim_to_data=True) - - channelGenericNoiseAdder.run(evt, station, det, type="rayleigh", amplitude=20 * units.mV) - - triggerSimulator.run(evt, station, det, number_concidences=2, threshold=100 * units.mV) - - if station.get_trigger('default_simple_threshold').has_triggered(): - - channelBandPassFilter.run(evt, station, det, passband=[80 * units.MHz, 500 * units.MHz], filter_type='butter', order=10) - - eventTypeIdentifier.run(evt, station, "forced", 'cosmic_ray') - - channelStopFilter.run(evt, station, det) - - channelBandPassFilter.run(evt, station, det, passband=[60 * units.MHz, 600 * units.MHz], filter_type='rectangular') - - channelSignalReconstructor.run(evt, station, det) - - hardwareResponseIncorporator.run(evt, station, det) - - correlationDirectionFitter.run(evt, station, det, n_index=1., channel_pairs=channel_pairs) - - voltageToEfieldConverter.run(evt, station, det, use_channels=used_channels_efield) - - electricFieldSignalReconstructor.run(evt, station, det) - voltageToAnalyticEfieldConverter.run(evt, station, det, use_channels=used_channels_efield, bandpass=[80 * units.MHz, 500 * units.MHz], use_MC_direction=False) - - channelResampler.run(evt, station, det, sampling_rate=1 * units.GHz) - - electricFieldResampler.run(evt, station, det, sampling_rate=1 * units.GHz) - - eventWriter.run(evt, det) - - -nevents = eventWriter.end() -print("Finished processing, {} events".format(nevents)) diff --git a/NuRadioReco/examples/cr_analysis/analyse_sim.py b/NuRadioReco/examples/cr_analysis/analyse_sim.py new file mode 100644 index 000000000..0e130e5fb --- /dev/null +++ b/NuRadioReco/examples/cr_analysis/analyse_sim.py @@ -0,0 +1,141 @@ +import numpy as np +import os +import glob +import helper_cr_eff as hcr +import json +import NuRadioReco.modules.io.eventReader as eventReader +from NuRadioReco.framework.parameters import stationParameters as stnp +from NuRadioReco.framework.parameters import showerParameters as shp +from NuRadioReco.utilities import units +import argparse +import logging +logger = logging.getLogger() +logger.setLevel(logging.INFO) + + +parser = argparse.ArgumentParser(description='Nurfile analyser') +parser.add_argument('--directory', type=str, nargs='?', + default='output_air_shower_reco/', + help='path were results from air shower analysis are stored') +parser.add_argument('--condition', type=str, nargs='?', default='', + help='string which should be in dict name') +parser.add_argument('--energy_binning', type=list, nargs='?', + default=[16.5, 20, 8], help='energy bins as log() with start, stop, number of bins (np.logspace)') +parser.add_argument('--zenith_binning', type=list, nargs='?', + default=[0, 100, 10], help='zenith bins in deg with start, stop, step (np.arange)') + +args = parser.parse_args() +energy_binning = args.energy_binning +zenith_binning = args.zenith_binning +distance_binning = args.distance_binning +number_of_sta_in_evt = args.number_of_sta_in_evt + +# the entries of this list are defined in the input argument energy_bins. +# [0] is the start value, [1] is the stop value, [2] is the number of samples generated +energy_bins = np.logspace(energy_binning[0], energy_binning[1], energy_binning[2]) +energy_bins_low = energy_bins[0:-1] +energy_bins_high = energy_bins[1:] +logger.info(f"Use energy bins {energy_bins} eV") + +# the entries of this list are defined in the input argument zenith_bins. +# [0] is the start value, [1] is the stop value, [2] is step size +zenith_bins = np.arange(zenith_binning[0], zenith_binning[1], zenith_binning[2]) * units.deg +zenith_bins_low = zenith_bins[0:-1] +zenith_bins_high = zenith_bins[1:] +logger.info(f"Use zenith bins {zenith_bins/units.deg} deg") + +nur_file_list = [] # get non corrupted input files with specified passband +i = 0 +for nur_file in glob.glob('{}*.nur'.format(args.directory)): + if os.path.isfile(nur_file) and str(args.condition) in nur_file: + i += 1 + if os.path.getsize(nur_file) > 0: + nur_file_list.append(nur_file) + +n_files = len(nur_file_list) + +energy = [] +zenith = [] +azimuth = [] +distance = [] +events = [] +trigger_status = [] # trigger status per event station and threshold +trigger_in_station = [] # name of trigger + +num = 0 + +evtReader = eventReader.eventReader() +evtReader.begin(filename=nur_file, read_detector=True) +det = evtReader.get_detector() +default_station = det.get_station_ids()[0] +for evt in evtReader.run(): # loop over all events, one event is one station + num += 1 + event_id = evt.get_id() + events.append(evt) + det_position = det.get_absolute_position(station_id=default_station) + sta = evt.get_station(station_id=default_station) + sim_station = sta.get_sim_station() + energy.append(sim_station.get_parameter(stnp.cr_energy)) + zenith.append(sim_station.get_parameter(stnp.zenith)) # get zenith for each station + azimuth.append(sim_station.get_parameter(stnp.azimuth)) + trigger_in_station.append(sta.get_triggers()) # get trigger for each station + trigger_status.append(sta.has_triggered()) + + for sho in evt.get_sim_showers(): + core = sho.get_parameter(shp.core) + distance.append(np.sqrt( + ((core[0] - det_position[0])**2) + + (core[1] - det_position[1])**2)) + +trigger_status = np.array(trigger_status) +zenith = np.array(zenith) +zenith_deg = zenith / units.deg # this is necessary to avoid mistakes due to decimals +distance = np.array(distance) +energy = np.array(energy) + + +# here we create empty array which will be filled in the following loop. The first axis contains all parameters in +# the energy bin and the second axis the zenthis bins + +triggered_trigger_e = np.zeros((len(energy_bins_low), len(zenith_bins_low))) +masked_events_e = np.zeros((len(energy_bins_low), len(zenith_bins_low))) +trigger_efficiency_e = np.zeros((len(energy_bins_low), len(zenith_bins_low))) + +for dim_0, energy_bin_low, energy_bin_high in zip(range(len(energy_bins_low)), energy_bins_low, energy_bins_high): + mask_e = (energy >= energy_bin_low) & (energy < energy_bin_high) # choose one energy bin + n_events_masked_e = np.sum(mask_e) # number of events in that energy bin + + for dim_1, zenith_bin_low, zenith_bin_high in zip(range(len(zenith_bins_low)), zenith_bins_low, zenith_bins_high): + # choose zenith bin + mask_z = (zenith_deg >= zenith_bin_low/units.deg) & (zenith_deg < zenith_bin_high/units.deg) + # trigger in in one energy bin and one zenith bin (ez) (values depend on the loop) + mask_ez = mask_e & mask_z + masked_trigger_status_ez = trigger_status[mask_ez] + n_events_masked_ez = np.sum(mask_ez) # number of events in that energy and zenith bin + + # number of triggered true trigger in energy and zentih bin + triggered_trigger_ez = np.sum(masked_trigger_status_ez, axis=0) + # fraction of events in this energy and zeniths bin that triggered true + trigger_efficiency_ez = triggered_trigger_ez / n_events_masked_ez + + +dic = {} +dic['default_station'] = default_station +dic['energy_bins_low'] = energy_bins_low +dic['energy_bins_high'] = energy_bins_high +dic['zenith_bins_low'] = zenith_bins_low +dic['zenith_bins_high'] = zenith_bins_high +dic['trigger_masked_events'] = np.nan_to_num(masked_events_e) +dic['triggered_trigger'] = np.nan_to_num(triggered_trigger_e) +dic['trigger_efficiency'] = np.nan_to_num(trigger_efficiency_e) + + +out_dir = 'results' +os.makedirs(out_dir, exist_ok=True) + +output_file = 'dict_air_shower_e{}_z{}.json'.format( + len(energy_bins_low), + len(zenith_bins_low)) + +with open(os.path.join(out_dir, output_file), 'w') as outfile: + json.dump(dic, outfile, cls=hcr.NumpyEncoder, indent=4) diff --git a/NuRadioReco/examples/cr_efficiency_analysis/helper_cr_eff.py b/NuRadioReco/examples/cr_analysis/helper_cr_eff.py similarity index 100% rename from NuRadioReco/examples/cr_efficiency_analysis/helper_cr_eff.py rename to NuRadioReco/examples/cr_analysis/helper_cr_eff.py diff --git a/NuRadioReco/examples/cr_efficiency_analysis/6_plot_aeff_cr_number.py b/NuRadioReco/examples/cr_analysis/plot_aeff_cr_num.py similarity index 81% rename from NuRadioReco/examples/cr_efficiency_analysis/6_plot_aeff_cr_number.py rename to NuRadioReco/examples/cr_analysis/plot_aeff_cr_num.py index e0612e640..9e7b44f03 100644 --- a/NuRadioReco/examples/cr_efficiency_analysis/6_plot_aeff_cr_number.py +++ b/NuRadioReco/examples/cr_analysis/plot_aeff_cr_num.py @@ -4,7 +4,7 @@ from NuRadioReco.utilities import units import NuRadioReco.utilities.cr_flux as hcr -filename = 'results/air_shower/dict_air_shower_pb_80_180_e7_z9_d2_4000.json' +filename = '' with open(filename, 'r') as fp: data = json.load(fp) @@ -18,17 +18,6 @@ trigger_effective_area_err = np.array(data['trigger_effective_area_err']) trigger_weight = np.array(data['triggered_trigger_weight']) -### choose only some energy and zenith bins for plotting -chosen_energy_bins = [0, -1] -chosen_zenith_bins = [0, 8] - -trigger_effective_area = trigger_effective_area[chosen_energy_bins[0]:chosen_energy_bins[1],chosen_zenith_bins[0]:chosen_zenith_bins[1]] -trigger_effective_area_err = trigger_effective_area_err[chosen_energy_bins[0]:chosen_energy_bins[1],chosen_zenith_bins[0]:chosen_zenith_bins[1]] -trigger_weight = trigger_weight[chosen_energy_bins[0]:chosen_energy_bins[1],chosen_zenith_bins[0]:chosen_zenith_bins[1]] -energy_bins_low = energy_bins_low[chosen_energy_bins[0]:chosen_energy_bins[1]] -energy_bins_high = energy_bins_high[chosen_energy_bins[0]:chosen_energy_bins[1]] -zenith_bins_low = zenith_bins_low[chosen_zenith_bins[0]:chosen_zenith_bins[1]] -zenith_bins_high = zenith_bins_high[chosen_zenith_bins[0]:chosen_zenith_bins[1]] steradian = [] weight_det = [] diff --git a/NuRadioReco/examples/cr_analysis/sim_cr_det_array.py b/NuRadioReco/examples/cr_analysis/sim_cr_det_array.py new file mode 100644 index 000000000..a6d4a39e1 --- /dev/null +++ b/NuRadioReco/examples/cr_analysis/sim_cr_det_array.py @@ -0,0 +1,102 @@ +import numpy as np +import NuRadioReco.examples.cr_efficiency_analysis.helper_cr_eff as hcr +from NuRadioReco.utilities import units +from NuRadioReco.detector.detector import Detector +from NuRadioReco.framework.trigger import RadiantAUXTrigger +import NuRadioReco.modules.io.coreas.readCoREASDetector +import NuRadioReco.modules.efieldToVoltageConverter +import NuRadioReco.modules.RNO_G.hardwareResponseIncorporator +import NuRadioReco.modules.channelGenericNoiseAdder +import NuRadioReco.modules.channelGalacticNoiseAdder_fast +import NuRadioReco.modules.trigger.radiant_aux_trigger +import NuRadioReco.modules.channelBandPassFilter +import NuRadioReco.modules.eventTypeIdentifier +import NuRadioReco.modules.channelResampler +import NuRadioReco.modules.electricFieldResampler +import NuRadioReco.modules.io.eventWriter +import logging +import argparse + +logger = logging.getLogger() +logger.setLevel(logging.WARNING) + +parser = argparse.ArgumentParser(description='Run air shower Reconstruction') + +parser.add_argument('--input_file', type=str, nargs='?', default='example_data/example_data.hdf5') +parser.add_argument('--det_file', type=str, nargs='?', default='../detector/RNO_G/RNO_season_2023.json') +args = parser.parse_args() +logger.info(f'Using detector file {args.det_file} on {args.input_file}') + +det = Detector(json_filename=args.det_file, antenna_by_depth=False) + +#core_positions = NuRadioReco.modules.io.coreas.readCoREASDetector.get_random_core_positions() + +# module to read the CoREAS file and convert it to NuRadioReco event for an array of detector stations. +readCoREASDetector = NuRadioReco.modules.io.coreas.readCoREASDetector.readCoREASDetector() +readCoREASDetector.begin(args.input_file, log_level=logging.WARNING) + +# module to set the event type, if cosmic ray, the refraction of the emission on the air ice interface is taken into account +eventTypeIdentifier = NuRadioReco.modules.eventTypeIdentifier.eventTypeIdentifier() + +# module to convolves the electric field with the antenna response +efieldToVoltageConverter = NuRadioReco.modules.efieldToVoltageConverter.efieldToVoltageConverter() +efieldToVoltageConverter.begin(debug=False) + +# module to add the detector response, e.g. amplifier, filter, etc. +hardwareResponseIncorporator = NuRadioReco.modules.RNO_G.hardwareResponseIncorporator.hardwareResponseIncorporator() + +# module to add thermal noise to the channels +channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() +channelGenericNoiseAdder.begin() + +# module to add galactic noise to the channels +channelGalacticNoiseAdder = NuRadioReco.modules.channelGalacticNoiseAdder_fast.channelGalacticNoiseAdder() +channelGalacticNoiseAdder.begin(skymodel='gsm2016', n_side=4, freq_range=np.array([0.07, 0.81])) + +# module to simulate the trigger +triggerSimulator = NuRadioReco.modules.trigger.radiant_aux_trigger.triggerSimulator() +triggerSimulator.begin() + +# module to filter the channels +channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() +channelBandPassFilter.begin() + +# adjusts sampling rate for electirc field +electricFieldResampler = NuRadioReco.modules.electricFieldResampler.electricFieldResampler() +electricFieldResampler.begin() + +# adjusts sampling rate for channels +channelResampler = NuRadioReco.modules.channelResampler.channelResampler() +channelResampler.begin() + +# module to write the events to a .nur file +eventWriter = NuRadioReco.modules.io.eventWriter.eventWriter() +eventWriter.begin(f'cr_reco_array.nur') + +for evt in readCoREASDetector.run(det): + for sta in evt.get_stations(): + eventTypeIdentifier.run(evt, sta, "forced", 'cosmic_ray') + + efieldToVoltageConverter.run(evt, sta, det) + + channelGenericNoiseAdder.run(evt, sta, det, amplitude=1.091242302378349e-05, min_freq=80*units.MHz, max_freq=800*units.MHz, type='rayleigh') + channelGalacticNoiseAdder.run(evt, sta, det, passband=[0.08, 0.8]) + + hardwareResponseIncorporator.run(evt, sta, det, sim_to_data=True) + + triggerSimulator.run(evt, sta, det, + threshold_sigma=10, + triggered_channels=[13, 16, 19], #surface channels + trigger_name=f'radiant_trigger_10sigma') + + channelResampler.run(evt, sta, det, sampling_rate=3.2) + + electricFieldResampler.run(evt, sta, det, sampling_rate=3.2) + + eventWriter.run(evt, det=det, mode={ + 'Channels': True, + 'ElectricFields': True, + 'SimChannels': False, + 'SimElectricFields': False + }) +nevents = eventWriter.end() \ No newline at end of file diff --git a/NuRadioReco/examples/cr_analysis/sim_cr_single_station.py b/NuRadioReco/examples/cr_analysis/sim_cr_single_station.py new file mode 100644 index 000000000..6e19d0b40 --- /dev/null +++ b/NuRadioReco/examples/cr_analysis/sim_cr_single_station.py @@ -0,0 +1,109 @@ +import numpy as np +import helper_cr_eff as hcr +from NuRadioReco.utilities import units +from NuRadioReco.detector.detector import Detector +import NuRadioReco.modules.io.coreas.readCoREASStation +import NuRadioReco.modules.efieldToVoltageConverter +import NuRadioReco.modules.RNO_G.hardwareResponseIncorporator +import NuRadioReco.modules.channelGenericNoiseAdder +import NuRadioReco.modules.channelGalacticNoiseAdder +import NuRadioReco.modules.trigger.highLowThreshold +import NuRadioReco.modules.channelBandPassFilter +import NuRadioReco.modules.eventTypeIdentifier +import NuRadioReco.modules.channelResampler +import NuRadioReco.modules.electricFieldResampler +import NuRadioReco.modules.io.eventWriter +import logging +import argparse + +logger = logging.getLogger() +logger.setLevel(logging.INFO) + + +parser = argparse.ArgumentParser(description='Run air shower Reconstruction') + +parser.add_argument('--detector_file', type=str, nargs='?', default='example_data/arianna_station_32.json', + help='choose detector with a single station for air shower simulation') +parser.add_argument('--input_file', type=str, nargs='?', + default='example_data/example_data.hdf5', help='hdf5 coreas file') + +args = parser.parse_args() + +logger.info(f"Use {args.detector_file} on file {args.input_file}") + +det = Detector(json_filename=args.detector_file) +station_id = det.get_station_ids()[0] + +# module to read the CoREAS file and convert it to NuRadioReco event, each observer is a new event with a different core position +readCoREASStation = NuRadioReco.modules.io.coreas.readCoREASStation.readCoREASStation() +readCoREASStation.begin([args.input_file], station_id, debug=False) + +# module to set the event type, if cosmic ray, the refraction of the emission on the air ice interface is taken into account +eventTypeIdentifier = NuRadioReco.modules.eventTypeIdentifier.eventTypeIdentifier() + +# module to convolves the electric field with the antenna response +efieldToVoltageConverter = NuRadioReco.modules.efieldToVoltageConverter.efieldToVoltageConverter() +efieldToVoltageConverter.begin(debug=False) + +# module to add the detector response, e.g. amplifier, filter, etc. +hardwareResponseIncorporator = NuRadioReco.modules.RNO_G.hardwareResponseIncorporator.hardwareResponseIncorporator() + +# module to add thermal noise to the channels +channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() +channelGenericNoiseAdder.begin() + +# module to add galactic noise to the channels +channelGalacticNoiseAdder = NuRadioReco.modules.channelGalacticNoiseAdder.channelGalacticNoiseAdder() +channelGalacticNoiseAdder.begin(n_side=4, interpolation_frequencies=np.arange(0.01, 0.81, 0.1)) + +# module to simulate the trigger +triggerSimulator = NuRadioReco.modules.trigger.highLowThreshold.triggerSimulator() +triggerSimulator.begin() + +# module to filter the channels +channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() +channelBandPassFilter.begin() + +# module adjust the sampling rate of the electric field +electricFieldResampler = NuRadioReco.modules.electricFieldResampler.electricFieldResampler() +electricFieldResampler.begin() + +# module adjust the sampling rate of the channels +channelResampler = NuRadioReco.modules.channelResampler.channelResampler() +channelResampler.begin() + +# module to write the event to a .nur file +eventWriter = NuRadioReco.modules.io.eventWriter.eventWriter() +eventWriter.begin('cr_single_station.nur') + +for evt in readCoREASStation.run(detector=det): + for sta in evt.get_stations(): + eventTypeIdentifier.run(evt, sta, "forced", 'cosmic_ray') + + efieldToVoltageConverter.run(evt, sta, det) + + channelGenericNoiseAdder.run(evt, sta, det, amplitude=1.091242302378349e-05, min_freq=80*units.MHz, max_freq=800*units.MHz, type='rayleigh') + + channelGalacticNoiseAdder.run(evt, sta, det) + + hardwareResponseIncorporator.run(evt, sta, det, sim_to_data=True) + + triggerSimulator.run(evt, sta, det, + threshold_high=5e-6, + threshold_low=-5e-6, + coinc_window=60, + number_concidences=2, + triggered_channels=[0, 1, 2, 3], + trigger_name='high_low') + + channelResampler.run(evt, sta, det, sampling_rate=3.2) + + electricFieldResampler.run(evt, sta, det, sampling_rate=3.2) + + eventWriter.run(evt, det=det, mode={ + 'Channels': True, + 'ElectricFields': True, + 'SimChannels': False, + 'SimElectricFields': False + }) +nevents = eventWriter.end() diff --git a/NuRadioReco/examples/cr_efficiency_analysis/1_create_config_file.py b/NuRadioReco/examples/cr_efficiency_analysis/1_create_config_file.py deleted file mode 100644 index da9417c31..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/1_create_config_file.py +++ /dev/null @@ -1,146 +0,0 @@ -import numpy as np -import helper_cr_eff as hcr -import json -import os -from NuRadioReco.utilities import units -import argparse - -parser = argparse.ArgumentParser() -parser.add_argument('--output_path', type=str, nargs='?', default=os.path.dirname(__file__), - help='Path to save output, most likely the path to the cr_efficiency_analysis directory') -parser.add_argument('--detector_file', type=str, nargs='?', - default='LPDA_Southpole.json', - help='file with one antenna at and the geographic location, change triggered channels accordingly') -parser.add_argument('--target_global_trigger_rate', type=float, nargs='?', default=100, - help='trigger rate for all channels in Hz') -parser.add_argument('--trigger_name', type=str, nargs='?', default='high_low', - help='name of the trigger, high_low, envelope or power_integration') -parser.add_argument('--default_station', type=int, nargs='?', default=101, - help='default station id') -parser.add_argument('--trace_samples', type=int, nargs='?', default=1024, - help='elements in the array of one trace') -parser.add_argument('--sampling_rate', type=int, nargs='?', default=1, - help='sampling rate in GHz') -parser.add_argument('--triggered_channels', type=np.ndarray, nargs='?', default=np.array([1]), - help='channel on which the trigger is applied') -parser.add_argument('--total_number_triggered_channels', type=int, nargs='?', default=4, - help='number ot channels that trigger.') -parser.add_argument('--number_coincidences', type=int, nargs='?', default=2, - help='number coincidences of true trigger within on station of the detector') -parser.add_argument('--coinc_window', type=int, nargs='?', default=80, - help='coincidence window within the number coincidence has to occur. In ns') -parser.add_argument('--int_window', type=float, nargs='?', default=10, - help='integration time window [ns] for power_integration trigger') -parser.add_argument('--passband_low', type=int, nargs='?', default=80, - help='lower bound of the passband used for the trigger in MHz') -parser.add_argument('--passband_high', type=int, nargs='?', default=180, - help='higher bound of the passband used for the trigger in MHz') -parser.add_argument('--order_trigger', type=int, nargs='?', default=10, - help='order of the filter used in the trigger') -parser.add_argument('--Tnoise', type=int, nargs='?', default=300, - help='Temperature of thermal noise in K') -parser.add_argument('--T_noise_min_freq', type=int, nargs='?', default=50, - help='min freq of thermal noise in MHz') -parser.add_argument('--T_noise_max_freq', type=int, nargs='?', default=800, - help='max freq of thermal noise in MHz') -parser.add_argument('--galactic_noise_n_side', type=int, nargs='?', default=4, - help='The n_side parameter of the healpix map. Has to be power of 2, basicly the resolution') -parser.add_argument('--galactic_noise_interpolation_frequencies_start', type=int, nargs='?', default=10, - help='start frequency the galactic noise is interpolated over in MHz') -parser.add_argument('--galactic_noise_interpolation_frequencies_stop', type=int, nargs='?', default=1100, - help='stop frequency the galactic noise is interpolated over in MHz') -parser.add_argument('--galactic_noise_interpolation_frequencies_step', type=int, nargs='?', default=100, - help='frequency steps the galactic noise is interpolated over in MHz') -parser.add_argument('--n_random_phase', type=int, nargs='?', default=10, - help='for computing time reasons one galactic noise amplitude is reused ' - 'n_random_phase times, each time a random phase is added') -parser.add_argument('--threshold_start', type=float, nargs='?', - help='value of the first tested threshold in Volt') -parser.add_argument('--threshold_step', type=float, nargs='?', - help='value of the threshold step in Volt') -parser.add_argument('--station_time', type=str, nargs='?', default='2021-01-01T00:00:00', - help='station time for calculation of galactic noise') -parser.add_argument('--station_time_random', type=bool, nargs='?', default=True, - help='choose if the station time should be random or not') -parser.add_argument('--hardware_response', type=bool, nargs='?', default=True, - help='choose if the hardware response (amp) should be True or False') -parser.add_argument('--iterations_per_job', type=int, nargs='?', default=200, - help='choose if the hardware response (amp) should be True or False') -parser.add_argument('--number_of_allowed_trigger', type=bool, nargs='?', default=3, - help='The number of iterations is calculated to yield a trigger rate') - -args = parser.parse_args() -target_global_trigger_rate = args.target_global_trigger_rate * units.Hz -passband_low = args.passband_low * units.megahertz -passband_high = args.passband_high * units.megahertz -passband_trigger = np.array([passband_low, passband_high]) -sampling_rate = args.sampling_rate * units.gigahertz -coinc_window = args.coinc_window * units.ns -int_window = args.int_window * units.ns -Tnoise = args.Tnoise * units.kelvin -T_noise_min_freq = args.T_noise_min_freq * units.megahertz -T_noise_max_freq = args.T_noise_max_freq * units.megahertz -galactic_noise_interpolation_frequencies_start = args.galactic_noise_interpolation_frequencies_start * units.MHz -galactic_noise_interpolation_frequencies_stop = args.galactic_noise_interpolation_frequencies_stop * units.MHz -galactic_noise_interpolation_frequencies_step = args.galactic_noise_interpolation_frequencies_step * units.MHz - -trace_length = args.trace_samples / sampling_rate -target_single_trigger_rate = hcr.get_single_channel_trigger_rate( - target_global_trigger_rate, args.total_number_triggered_channels, args.number_coincidences, coinc_window) -n_iteration_for_one_allowed_trigger = (trace_length * target_single_trigger_rate) ** -1 - -n_iterations = int(n_iteration_for_one_allowed_trigger * args.number_of_allowed_trigger / args.n_random_phase) -resolution = (n_iteration_for_one_allowed_trigger * args.number_of_allowed_trigger * trace_length) ** -1 - -number_of_jobs = n_iterations / args.iterations_per_job -number_of_jobs = int(np.ceil(number_of_jobs)) -Vrms_thermal_noise = hcr.calculate_thermal_noise_Vrms(Tnoise, T_noise_max_freq, T_noise_min_freq) - -if args.threshold_start is None: - if args.trigger_name == "power_integration": - raise Exception('Please set threshold start value manually for the power integration trigger') - else: - if args.hardware_response: - threshold_start = 1e3 * Vrms_thermal_noise - elif not args.hardware_response: - threshold_start = 1.8 * Vrms_thermal_noise -else: - threshold_start = args.threshold_start * units.volt - -if args.threshold_step is None: - if args.trigger_name == "power_integration": - raise Exception('Please set threshold step value manually for the power integration trigger') - else: - if args.hardware_response: - threshold_step = 1e-3 * units.volt - elif not args.hardware_response: - threshold_step = 1e-6 * units.volt -else: - threshold_step = args.threshold_step * units.volt - -dic = {'T_noise': Tnoise, 'Vrms_thermal_noise': Vrms_thermal_noise, 'n_iterations_total': n_iterations, - 'number_of_allowed_trigger': args.number_of_allowed_trigger, 'iterations_per_job': args.iterations_per_job, - 'number_of_jobs': number_of_jobs, 'target_single_trigger_rate': target_single_trigger_rate, - 'target_global_trigger_rate': target_global_trigger_rate, 'resolution': resolution, - 'trigger_name': args.trigger_name, 'passband_trigger': passband_trigger, - 'total_number_triggered_channels': args.total_number_triggered_channels, - 'number_coincidences': args.number_coincidences, 'triggered_channels': args.triggered_channels, - 'coinc_window': coinc_window, 'int_window': int_window, 'order_trigger': args.order_trigger, 'detector_file': args.detector_file, - 'default_station': args.default_station, 'trace_samples': args.trace_samples, 'sampling_rate': sampling_rate, - 'trace_length': trace_length, 'T_noise_min_freq': T_noise_min_freq, 'T_noise_max_freq': T_noise_max_freq, - 'galactic_noise_n_side': args.galactic_noise_n_side, - 'galactic_noise_interpolation_frequencies_start': galactic_noise_interpolation_frequencies_start, - 'galactic_noise_interpolation_frequencies_stop': galactic_noise_interpolation_frequencies_stop, - 'galactic_noise_interpolation_frequencies_step': galactic_noise_interpolation_frequencies_step, - 'station_time': args.station_time, 'station_time_random': args.station_time_random, - 'hardware_response': args.hardware_response, 'n_random_phase': args.n_random_phase, - 'threshold_start': threshold_start, 'threshold_step': threshold_step} - -os.makedirs(os.path.join(args.output_path, 'config/ntr'), exist_ok=True) - -output_file = f'config/ntr/config_{args.trigger_name}_trigger_rate_{target_global_trigger_rate/units.Hz:.0f}Hz_coinc_{args.number_coincidences}of{args.total_number_triggered_channels}.json' - -abs_path_output_file = os.path.normpath(os.path.join(args.output_path, output_file)) - -with open(abs_path_output_file, 'w') as outfile: - json.dump(dic, outfile, cls=hcr.NumpyEncoder, indent=4, sort_keys=True) diff --git a/NuRadioReco/examples/cr_efficiency_analysis/2_II_merge_output_calculate_trigger_rate.py b/NuRadioReco/examples/cr_efficiency_analysis/2_II_merge_output_calculate_trigger_rate.py deleted file mode 100644 index 1f6e81d08..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/2_II_merge_output_calculate_trigger_rate.py +++ /dev/null @@ -1,84 +0,0 @@ -import numpy as np -import helper_cr_eff as hcr -import matplotlib.pyplot as plt -import os -from NuRadioReco.utilities import units -import argparse -import json -import glob -import logging -logger = logging.getLogger() -logger.setLevel(logging.INFO) - -''' -Use the script calculate trigger_rate_for_threshold.py to create the files. -This script reads all files in output_treshold_calculations and writes one new dict which then contains the -information of all files. This is necessary if you want to split the total number of iterations in different job. -''' - -parser = argparse.ArgumentParser(description='Noise Trigger Rate') -parser.add_argument('--directory', type=str, nargs='?', - default='output_threshold_calculation/', help='directory with output files') -parser.add_argument('--condition', type=str, nargs='?', - default='envelope_trigger_0Hz_3of4', help='string which should be in dict name') -args = parser.parse_args() - -logger.info(f"Checking {args.directory} for condition: {args.condition}") - -file_list = [] -# get non corrupted files from threshold calculations with specified passband -for file in glob.glob(args.directory + '*' + args.condition + '*.json'): - if os.path.isfile(file) and os.path.getsize(file) > 0: - file_list.append(file) - - -n_files = len(file_list) - -logger.info(f"Using files {file_list}") - -# open one file to check the number of tested thresholds -with open(file_list[0], 'r') as fp: - cfg = json.load(fp) - -thresholds = cfg['thresholds'] - -trigger_efficiency = np.zeros((n_files, len(thresholds))) -trigger_rate = np.zeros((n_files, len(thresholds))) - -for i_file, filename in enumerate(file_list): - with open(filename, 'r') as fp: - data = json.load(fp) - trigger_efficiency[i_file] = data['efficiency'] - trigger_rate[i_file] = data['trigger_rate'] - -trigger_efficiency_all = np.sum(trigger_efficiency, axis=0) / n_files -trigger_rate_all = np.sum(trigger_rate, axis=0) / n_files -iterations = cfg['n_iterations_total'] * cfg['n_random_phase'] * n_files -estimated_global_rate = hcr.get_global_trigger_rate(trigger_rate_all, cfg['total_number_triggered_channels'], - cfg['number_coincidences'], cfg['coinc_window']) - -dic = {} -dic = cfg.copy() -dic['iteration'] = iterations -dic['efficiency'] = trigger_efficiency_all -dic['trigger_rate'] = trigger_rate_all -dic['estimated_global_trigger_rate'] = estimated_global_rate - -final_threshold = thresholds[np.argmin(trigger_rate_all - cfg['target_single_trigger_rate'])] - -dic['final_threshold'] = final_threshold - -os.makedirs('config/air_shower', exist_ok=True) - -if cfg['hardware_response']: - output_file = 'config/air_shower/final_config_{}_trigger_{:.0f}Hz_{}of{}_{:.2f}mV.json'.format( - cfg['trigger_name'], cfg['target_global_trigger_rate']/units.Hz, - cfg['number_coincidences'], cfg['total_number_triggered_channels'], - final_threshold/units.millivolt) -else: - output_file = 'config/air_shower/final_config_{}_trigger_{:.0f}Hz_{}of{}_{:.2e}V.json'.format( - cfg['trigger_name'], cfg['target_global_trigger_rate']/units.Hz, - cfg['number_coincidences'], cfg['total_number_triggered_channels'], final_threshold) - -with open(output_file, 'w') as outfile: - json.dump(dic, outfile, cls=hcr.NumpyEncoder, indent=4, sort_keys=True) diff --git a/NuRadioReco/examples/cr_efficiency_analysis/2_I_calculate_trigger_rate_for_threshold.py b/NuRadioReco/examples/cr_efficiency_analysis/2_I_calculate_trigger_rate_for_threshold.py deleted file mode 100644 index 73f6c8770..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/2_I_calculate_trigger_rate_for_threshold.py +++ /dev/null @@ -1,202 +0,0 @@ -import numpy as np -import helper_cr_eff as hcr -import json -import os -import time -import NuRadioReco.modules.channelGenericNoiseAdder -import NuRadioReco.modules.channelGalacticNoiseAdder -import NuRadioReco.modules.trigger.envelopeTrigger as envelopeTrigger -import NuRadioReco.modules.trigger.highLowThreshold as highLowThreshold -import NuRadioReco.modules.trigger.powerIntegration as powerIntegration -import NuRadioReco.modules.RNO_G.hardwareResponseIncorporator -import NuRadioReco.modules.channelBandPassFilter -import NuRadioReco.modules.eventTypeIdentifier -import NuRadioReco.utilities.fft -from NuRadioReco.detector.generic_detector import GenericDetector -from NuRadioReco.utilities import units -import argparse -import logging - -logger = logging.getLogger() -logger.setLevel(logging.INFO) - -''' -This script calculates the trigger rate for different thresholds over a given number of iterations. -The resolution is calculated in the create_config_file. - -For the galactic noise, the sky maps from PyGDSM are used. You can install it with -pip install git+https://github.com/telegraphic/pygdsm . - -The sampling rate has a huge influence on the threshold, because the trace has more time to exceed the threshold -for a sampling rate of 1GHz, 1955034 iterations yields a resolution of 0.5 Hz -Due to computational efficiency (galactic noise adder is slow), one amplitude is reused with 10 random phases - -For one config file I used something on the cluster like: -for number in $(seq 0 1 110); do echo $number; qsub Cluster_ntr_2.sh --number $number; sleep 0.2; done; - -for several config files I used: -FILES=output_threshold_calculation/* -for file in $FILES -do -echo "Processing $file" - for number in $(seq 0 1 110) - do echo $number - qsub calculate_threshold_multi_job --config_file $file --number $number - sleep 0.2 - done -done -''' - -parser = argparse.ArgumentParser(description='Noise Trigger Rate') -parser.add_argument('--config_file', type=str, nargs='?', - help='input filename from which the calculation starts.') -parser.add_argument('--number', type=int, nargs='?', default=1, - help='Assigns a number to the scrip. ' - 'Check in config file for number_of_jobs ' - 'and sumbit the job with a for loop over the range of number_of_jobs') -parser.add_argument('--n_thresholds', type=int, nargs='?', default=15, - help='number of thresholds that will be tested.') -parser.add_argument('--output_path', type=str, nargs='?', default=os.path.dirname(__file__), - help='Path to save output, most likely the path to the cr_efficiency_analysis directory') - -args = parser.parse_args() - -with open(args.config_file, 'r') as fp: - cfg = json.load(fp) - -n_iterations_total = cfg['n_iterations_total'] -trigger_thresholds = (np.arange(cfg['threshold_start'], - cfg['threshold_start'] + (args.n_thresholds * cfg['threshold_step']), - cfg['threshold_step'])) - -logger.info("Processing trigger thresholds {}".format(trigger_thresholds)) - -det = GenericDetector(json_filename=cfg['detector_file'], default_station=cfg['default_station']) -station_ids = det.get_station_ids() -channel_ids = det.get_channel_ids(station_ids[0]) - -event, station, channel = hcr.create_empty_event( - det, cfg['trace_samples'], cfg['station_time'], cfg['station_time_random'], cfg['sampling_rate']) - -eventTypeIdentifier = NuRadioReco.modules.eventTypeIdentifier.eventTypeIdentifier() - -channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() -channelGenericNoiseAdder.begin() - -channelGalacticNoiseAdder = NuRadioReco.modules.channelGalacticNoiseAdder.channelGalacticNoiseAdder() -channelGalacticNoiseAdder.begin( - n_side=cfg['galactic_noise_n_side'], - interpolation_frequencies=np.arange(cfg['galactic_noise_interpolation_frequencies_start'], - cfg['galactic_noise_interpolation_frequencies_stop'], - cfg['galactic_noise_interpolation_frequencies_step'])) - -hardwareResponseIncorporator = NuRadioReco.modules.RNO_G.hardwareResponseIncorporator.hardwareResponseIncorporator() - -channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() -channelBandPassFilter.begin() - -t = time.time() # absolute time of system -sampling_rate = station.get_channel(channel_ids[0]).get_sampling_rate() -dt = 1. / sampling_rate - -time = channel.get_times() -channel_trace_start_time = time[0] -channel_trace_final_time = time[len(time) - 1] -channel_trace_time_interval = channel_trace_final_time - channel_trace_start_time - -trigger_status = [] -trigger_rate = [] -trigger_efficiency = [] - -for n_it in range(n_iterations_total): - station = event.get_station(cfg['default_station']) - if cfg['station_time_random']: - station = hcr.set_random_station_time(station, cfg['station_time']) - - eventTypeIdentifier.run(event, station, "forced", 'cosmic_ray') - - channel = hcr.create_empty_channel_trace(station, cfg['trace_samples'], sampling_rate) - - channelGenericNoiseAdder.run(event, station, det, - amplitude=cfg['Vrms_thermal_noise'], - min_freq=cfg['T_noise_min_freq'], - max_freq=cfg['T_noise_max_freq'], type='rayleigh', bandwidth=None) - - channelGalacticNoiseAdder.run(event, station, det) - - if cfg['hardware_response']: - hardwareResponseIncorporator.run(event, station, det, sim_to_data=True) - - for i_phase in range(cfg['n_random_phase']): - channel = hcr.add_random_phase(station, sampling_rate) - channelBandPassFilter.run(event, station, det, passband=cfg['passband_trigger'], - filter_type='butter', order=cfg['order_trigger']) - - trace = station.get_channel(station.get_channel_ids()[0]).get_trace() - trigger_status_all_thresholds = [] - for threshold in trigger_thresholds: - if cfg['trigger_name'] == 'high_low': - triggered_samples = highLowThreshold.get_high_low_triggers(trace, threshold, -threshold, - cfg['coinc_window'], dt=dt) - if True in triggered_samples: - has_triggered = bool(1) - else: - has_triggered = bool(0) - - if cfg['trigger_name'] == 'envelope': - triggered_samples = envelopeTrigger.get_envelope_triggers(trace, threshold) - if True in triggered_samples: - has_triggered = bool(1) - else: - has_triggered = bool(0) - - if cfg['trigger_name'] == 'power_integration': - trace_1 = station.get_channel(1).get_trace() - triggered_samples, int_power = powerIntegration.get_power_int_triggers(trace, threshold, cfg['int_window'], dt=dt, full_output=True) - if True in triggered_samples: - has_triggered = bool(1) - else: - has_triggered = bool(0) - print('threshold', threshold) - print('int_power', np.max(int_power)) - print(has_triggered) - - trigger_status_all_thresholds.append(has_triggered) - trigger_status.append(trigger_status_all_thresholds) - -trigger_status = np.array(trigger_status) -triggered_trigger = np.sum(trigger_status, axis=0) -trigger_efficiency = triggered_trigger / len(trigger_status) -trigger_rate = (1 / channel_trace_time_interval) * trigger_efficiency - -logger.info("Triggered true per trigger thresholds {}/{}".format(triggered_trigger, len(trigger_status))) - -dic = {'thresholds': trigger_thresholds, 'efficiency': trigger_efficiency, 'trigger_rate': trigger_rate, - 'T_noise': cfg['T_noise'], 'Vrms_thermal_noise': cfg['Vrms_thermal_noise'], - 'n_iterations_total': n_iterations_total, 'iterations_per_job': cfg['iterations_per_job'], - 'number_of_jobs': cfg['number_of_jobs'], 'target_single_trigger_rate': cfg['target_single_trigger_rate'], - 'target_global_trigger_rate': cfg['target_global_trigger_rate'], 'resolution': cfg['resolution'], - 'trigger_name': cfg['trigger_name'], 'passband_trigger': cfg['passband_trigger'], - 'total_number_triggered_channels': cfg['total_number_triggered_channels'], - 'number_coincidences': cfg['number_coincidences'], 'triggered_channels': cfg['triggered_channels'], - 'coinc_window': cfg['coinc_window'], 'order_trigger': cfg['order_trigger'], - 'detector_file': cfg['detector_file'], 'default_station': cfg['default_station'], - 'trace_samples': cfg['trace_samples'], 'sampling_rate': cfg['sampling_rate'], - 'trace_length': cfg['trace_length'], 'T_noise_min_freq': cfg['T_noise_min_freq'], - 'T_noise_max_freq': cfg['T_noise_max_freq'], 'galactic_noise_n_side': cfg['galactic_noise_n_side'], - 'galactic_noise_interpolation_frequencies_start': cfg['galactic_noise_interpolation_frequencies_start'], - 'galactic_noise_interpolation_frequencies_stop': cfg['galactic_noise_interpolation_frequencies_stop'], - 'galactic_noise_interpolation_frequencies_step': cfg['galactic_noise_interpolation_frequencies_step'], - 'station_time': cfg['station_time'], 'station_time_random': cfg['station_time_random'], - 'hardware_response': cfg['hardware_response'], 'n_random_phase': cfg['n_random_phase'], - 'threshold_start': cfg['threshold_start'], 'threshold_step': cfg['threshold_step']} - -os.makedirs(os.path.join(args.output_path, 'output_threshold_calculation'), exist_ok=True) - -output_file = 'output_threshold_calculation/{}_trigger_{:.0f}Hz_{}of{}_i{}_{}.json'.format( - cfg['trigger_name'], cfg['target_global_trigger_rate'] / units.Hz, - cfg['number_coincidences'], cfg['total_number_triggered_channels'], len(trigger_status), args.number) - -abs_path_output_file = os.path.normpath(os.path.join(args.output_path, output_file)) -with open(abs_path_output_file, 'w') as outfile: - json.dump(dic, outfile, cls=hcr.NumpyEncoder, indent=4, sort_keys=True) diff --git a/NuRadioReco/examples/cr_efficiency_analysis/2_calculate_threshold_for_trigger_rate.py b/NuRadioReco/examples/cr_efficiency_analysis/2_calculate_threshold_for_trigger_rate.py deleted file mode 100644 index 2ef367dc0..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/2_calculate_threshold_for_trigger_rate.py +++ /dev/null @@ -1,237 +0,0 @@ -import numpy as np -import helper_cr_eff as hcr -import time -import json -import os -import NuRadioReco.modules.channelGenericNoiseAdder -import NuRadioReco.modules.channelGalacticNoiseAdder -import NuRadioReco.modules.channelBandPassFilter -import NuRadioReco.modules.trigger.envelopeTrigger as envelopeTrigger -import NuRadioReco.modules.trigger.highLowThreshold as highLowThreshold -import NuRadioReco.modules.trigger.powerIntegration as powerIntegration -import NuRadioReco.modules.RNO_G.hardwareResponseIncorporator -import NuRadioReco.modules.eventTypeIdentifier -from NuRadioReco.detector.generic_detector import GenericDetector -from NuRadioReco.utilities import units -import argparse -import logging - -logger = logging.getLogger() -logger.setLevel(logging.INFO) - -''' -This script calculates the trigger threshold for a given global trigger rate. -The trigger rate for a single antenna is calculated in the crate_config_file.py -Afterwards the threshold is increased incrementally until the target trigger rate is achieved. -This script is slower than 2_I_calculate_trigger_rate_for_threshold, where the trigger rate -is calculated for a given array of thresholds. - -For the galactic noise, the sky maps from PyGDSM are used. You can install it with -pip install git+https://github.com/telegraphic/pygdsm . - -The sampling rate has a huge influence on the threshold, because the trace has more time to exceed the threshold -for a sampling rate of 1GHz, 1955034 iterations yields a resolution of 0.5 Hz -Due to computational efficiency (galactic noise adder is slow), one amplitude is reused with 10 random phases -''' - -parser = argparse.ArgumentParser(description='Noise Trigger Rate') -parser.add_argument('--config_file', type=str, nargs='?', - help='input filename from which the calculation starts.') -parser.add_argument('--output_path', type=str, nargs='?', default=os.path.dirname(__file__), - help='Path to save output, most likely the path to the cr_efficiency_analysis directory') - -args = parser.parse_args() - -with open(args.config_file, 'r') as fp: - cfg = json.load(fp) - -det = GenericDetector(json_filename=cfg['detector_file'], default_station=cfg['default_station']) -station_ids = det.get_station_ids() -channel_ids = det.get_channel_ids(station_ids[0]) - -logger.info("Apply {} trigger with passband {} " - .format(cfg['trigger_name'], np.array(cfg['passband_trigger']) / units.MHz)) - -event, station, channel = hcr.create_empty_event( - det, cfg['trace_samples'], cfg['station_time'], cfg['station_time_random'], cfg['sampling_rate']) - -eventTypeIdentifier = NuRadioReco.modules.eventTypeIdentifier.eventTypeIdentifier() -channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() -channelGenericNoiseAdder.begin() -channelGalacticNoiseAdder = NuRadioReco.modules.channelGalacticNoiseAdder.channelGalacticNoiseAdder() -channelGalacticNoiseAdder.begin( - n_side=cfg['galactic_noise_n_side'], - interpolation_frequencies=np.arange(cfg['galactic_noise_interpolation_frequencies_start'], - cfg['galactic_noise_interpolation_frequencies_stop'], - cfg['galactic_noise_interpolation_frequencies_step'])) -hardwareResponseIncorporator = NuRadioReco.modules.RNO_G.hardwareResponseIncorporator.hardwareResponseIncorporator() -channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() -channelBandPassFilter.begin() - -t = time.time() # absolute time of system -sampling_rate = station.get_channel(channel_ids[0]).get_sampling_rate() -dt = 1. / sampling_rate - -time = station.get_channel(channel_ids[0]).get_times() -channel_trace_start_time = time[0] -channel_trace_final_time = time[len(time) - 1] -channel_trace_time_interval = channel_trace_final_time - channel_trace_start_time - -trigger_status = [] -triggered_trigger = [] -trigger_rate = [] -trigger_efficiency = [] -thresholds = [] -iterations = [] -channel_rms = [] -channel_sigma = [] - -n_thres = 0 -sum_trigger = cfg['number_of_allowed_trigger'] + 1 -while sum_trigger > cfg['number_of_allowed_trigger']: - # with each iteration the threshold increases one step - threshold = cfg['threshold_start'] + (n_thres * cfg['threshold_step']) - thresholds.append(threshold) - logger.info("Processing threshold {:.2e} V".format(threshold)) - trigger_status_per_all_it = [] - - for n_it in range(cfg['n_iterations_total']): - station = event.get_station(cfg['default_station']) - eventTypeIdentifier.run(event, station, "forced", 'cosmic_ray') - - # here an empty channel trace is created - channel = hcr.create_empty_channel_trace(station, cfg['trace_samples'], sampling_rate) - - # thermal and galactic noise is added - channelGenericNoiseAdder.run(event, station, det, amplitude=cfg['Vrms_thermal_noise'], - min_freq=cfg['T_noise_min_freq'], max_freq=cfg['T_noise_max_freq'], - type='rayleigh') - channelGalacticNoiseAdder.run(event, station, det) - - # includes the amplifier response, if set true at the beginning - if cfg['hardware_response']: - hardwareResponseIncorporator.run(event, station, det, sim_to_data=True) - - channelBandPassFilter.run(event, station, det, passband=cfg['passband_trigger'], - filter_type='butter', order=cfg['order_trigger']) - - # This loop changes the phase of a trace with rand_phase, this is because the GalacticNoiseAdder - # needs some time and one amplitude is good enough for several traces. - # The current number of iteration can be calculated with i_phase + n_it*n_random_phase - for i_phase in range(cfg['n_random_phase']): - trigger_status_one_it = [] - channel = hcr.add_random_phase(station, sampling_rate) - - trace = station.get_channel(station.get_channel_ids()[0]).get_trace() - - if cfg['trigger_name'] == 'high_low': - triggered_samples = highLowThreshold.get_high_low_triggers(trace, threshold, -threshold, - cfg['coinc_window'], dt=dt) - if True in triggered_samples: - has_triggered = bool(1) - else: - has_triggered = bool(0) - - if cfg['trigger_name'] == 'envelope': - triggered_samples = envelopeTrigger.get_envelope_triggers(trace, threshold) - if True in triggered_samples: - has_triggered = bool(1) - else: - has_triggered = bool(0) - - if cfg['trigger_name'] == 'power_integration': - triggered_samples = powerIntegration.get_power_int_triggers(trace, threshold, cfg['int_window'], dt=dt, full_output=False) - if True in triggered_samples: - has_triggered = bool(1) - else: - has_triggered = bool(0) - - # trigger status for all iteration - trigger_status_per_all_it.append(has_triggered) - sum_trigger = np.sum(trigger_status_per_all_it) - - # here it is checked, how many of the triggers in n_iteration are triggered true. - # If it is more than number of allowed trigger, the threshold is increased with n_thres. - if np.sum(trigger_status_per_all_it) > cfg['number_of_allowed_trigger']: - number_of_trigger = np.sum(trigger_status_per_all_it) - trigger_efficiency_per_tt = np.sum(trigger_status_per_all_it) / len(trigger_status_per_all_it) - trigger_rate_per_tt = (1 / channel_trace_time_interval) * trigger_efficiency_per_tt - - trigger_rate.append(trigger_rate_per_tt) - trigger_efficiency.append(trigger_efficiency_per_tt) - estimated_global_rate = hcr.get_global_trigger_rate(trigger_rate_per_tt, cfg['total_number_triggered_channels'], - cfg['number_coincidences'], cfg['coinc_window']) - logger.info("current trigger rate {:.2e} Hz at threshold {:.2e} V" - .format(trigger_rate_per_tt / units.Hz, threshold)) - logger.info("target trigger rate single:{:.2e} Hz, global:{:.2e} Hz " - .format(cfg['target_single_trigger_rate'] / units.Hz, cfg['target_global_trigger_rate'] / units.Hz)) - logger.info("continue".format(cfg['trigger_name'])) - n_thres += 1 - - elif n_it == (cfg['n_iterations_total'] - 1): - number_of_trigger = np.sum(trigger_status_per_all_it) - trigger_efficiency_per_tt = np.sum(trigger_status_per_all_it) / len(trigger_status_per_all_it) - trigger_rate_per_tt = (1 / channel_trace_time_interval) * trigger_efficiency_per_tt - - trigger_rate.append(trigger_rate_per_tt) - trigger_efficiency.append(trigger_efficiency_per_tt) - estimated_global_rate = hcr.get_global_trigger_rate(trigger_rate_per_tt, cfg['total_number_triggered_channels'], - cfg['number_coincidences'], cfg['coinc_window']) - logger.info("checked {} thresholds".format(n_thres)) - logger.info("current trigger rate {:.2e} Hz at threshold {:.2e} V" - .format(trigger_rate_per_tt / units.Hz, threshold)) - logger.info("resolution for single trigger rate {:.2e} Hz" - .format(cfg['resolution'] / units.Hz)) - logger.info("estimated global trigger rate {:.2e} Hz" - .format(estimated_global_rate / units.Hz)) - logger.info("target trigger rate single:{:.2e} Hz, global:{:.2e} Hz" - .format(cfg['target_single_trigger_rate'] / units.Hz, cfg['target_global_trigger_rate'] / units.Hz)) - - dic = {'thresholds': thresholds, - 'efficiency': trigger_efficiency, - 'trigger_rate': trigger_rate, - 'estimated_global_trigger_rate': estimated_global_rate, - 'final_threshold': thresholds[-1], - 'T_noise': cfg['T_noise'], - 'Vrms_thermal_noise': cfg['Vrms_thermal_noise'], - 'n_iterations_total': cfg['n_iterations_total'], - 'iterations_per_job': cfg['iterations_per_job'], - 'number_of_jobs': cfg['number_of_jobs'], - 'target_single_trigger_rate': cfg['target_single_trigger_rate'], - 'target_global_trigger_rate': cfg['target_global_trigger_rate'], - 'resolution': cfg['resolution'], - 'trigger_name': cfg['trigger_name'], - 'passband_trigger': cfg['passband_trigger'], - 'total_number_triggered_channels': cfg['total_number_triggered_channels'], - 'number_coincidences': cfg['number_coincidences'], - 'triggered_channels': cfg['triggered_channels'], - 'coinc_window': cfg['coinc_window'], - 'order_trigger': cfg['order_trigger'], - 'detector_file': cfg['detector_file'], - 'default_station': cfg['default_station'], - 'trace_samples': cfg['trace_samples'], - 'sampling_rate': cfg['sampling_rate'], - 'trace_length': cfg['trace_length'], - 'T_noise_min_freq': cfg['T_noise_min_freq'], - 'T_noise_max_freq ': cfg['T_noise_max_freq'], - 'galactic_noise_n_side': cfg['galactic_noise_n_side'], - 'galactic_noise_interpolation_frequencies_start': cfg['galactic_noise_interpolation_frequencies_start'], - 'galactic_noise_interpolation_frequencies_stop': cfg['galactic_noise_interpolation_frequencies_stop'], - 'galactic_noise_interpolation_frequencies_step': cfg['galactic_noise_interpolation_frequencies_step'], - 'station_time': cfg['station_time'], - 'station_time_random': cfg['station_time_random'], - 'hardware_response': cfg['hardware_response'], - 'n_random_phase': cfg['n_random_phase'], - 'threshold_start': cfg['threshold_start'], - 'threshold_step': cfg['threshold_step'] - } - - os.makedirs(os.path.join(args.output_path, 'config/air_shower'), exist_ok=True) - - output_file = 'config/air_shower/final_config_{}_trigger_{:.2e}_{}of{}_{:.0f}Hz.json'.format( - cfg['trigger_name'], thresholds[-1], cfg['target_global_trigger_rate']/ units.Hz, - cfg['number_coincidences'], cfg['total_number_triggered_channels']) - abs_path_output_file = os.path.normpath(os.path.join(args.output_path, output_file)) - - with open(abs_path_output_file, 'w') as outfile: - json.dump(dic, outfile, cls=hcr.NumpyEncoder, indent=4, sort_keys=True) diff --git a/NuRadioReco/examples/cr_efficiency_analysis/3_plot_trigger_rate_vs_threshold.py b/NuRadioReco/examples/cr_efficiency_analysis/3_plot_trigger_rate_vs_threshold.py deleted file mode 100644 index 7c5fbccfe..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/3_plot_trigger_rate_vs_threshold.py +++ /dev/null @@ -1,92 +0,0 @@ -import numpy as np -import helper_cr_eff as hcr -import matplotlib.pyplot as plt -import os -from NuRadioReco.utilities import units -import argparse -import json - -'''This scripts plots the trigger rate for different thresholds -as given by the final config files for the air shower calculations. -''' - -parser = argparse.ArgumentParser(description='Noise Trigger Rate') -parser.add_argument('--config_directory', type=str, nargs='?', default='config/air_shower_envelope_trigger_1Hz_2of3_22.46mV/', - help='directory with final config files') -parser.add_argument('--save_plot', type=bool, nargs='?', default=True, help='decide if you want to save the plot') - -args = parser.parse_args() - -for config_file in os.listdir(args.config_directory): - with open(os.path.join(args.config_directory, config_file), 'r') as fp: - cfg = json.load(fp) - from scipy.interpolate import interp1d - - thresholds = np.array(cfg['thresholds']) - - # choose between a single antenna and global trigger rate - trigger_rate = np.array(cfg['estimated_global_trigger_rate']) - - - if cfg['hardware_response']: - x = thresholds / units.millivolt - y = trigger_rate / units.Hz - f1 = interp1d(x, y, 'cubic') - - xnew = np.linspace(thresholds[0] / units.millivolt, - thresholds[-1] / units.millivolt) - xlabel = 'Threshold [mV]' - target_trigger_rate = cfg['target_single_trigger_rate'] / units.Hz - - else: - x = thresholds / units.microvolt - y = trigger_rate / units.Hz - f1 = interp1d(x, y, 'cubic') - - xnew = np.linspace(thresholds[0] / units.microvolt, - thresholds[-1] / units.microvolt) - xlabel = r'Threshold [$\mu$V]' - target_trigger_rate = cfg['target_single_trigger_rate'] / units.Hz - - ynew = f1(xnew) - thresh = xnew[np.argmin(ynew - target_trigger_rate)] - - - plt.plot(x, y, marker='x', label='Noise trigger rate', - linestyle='none') - plt.plot(xnew, f1(xnew), '--', label='interp1d f({:.2e}) = {:.2e} Hz'.format(thresh, f1(thresh))) - plt.title('{} trigger, {}/{}, trigger rate {} Hz'.format(cfg['trigger_name'], - cfg['number_coincidences'], - cfg['total_number_triggered_channels'], - cfg['target_global_trigger_rate'] / units.Hz)) - plt.xlabel(xlabel, fontsize=18) - plt.hlines(target_trigger_rate, x[0], x[-1], label='Target trigger rate one antenna {:.2e} Hz'.format( - cfg['target_single_trigger_rate'] / units.Hz)) - plt.ylabel('Trigger rate [Hz]', fontsize=18) - plt.tick_params(axis='x', labelsize=16) - plt.tick_params(axis='y', labelsize=16) - plt.yscale('log') - plt.ylim(target_trigger_rate * 1e-2) - plt.legend() - plt.tight_layout() - - if args.save_plot: - os.makedirs('results/plots', exist_ok=True) - if cfg['hardware_response']: - plt.savefig( - 'results/plots/fig_ntr_{}_{:.0f}Hz_{}of{}_threshold_{:.2f}mV.png'.format( - cfg['trigger_name'], - cfg['target_global_trigger_rate'] / units.Hz, - cfg['number_coincidences'], - cfg['total_number_triggered_channels'], - cfg['final_threshold'] / units.millivolt)) - else: - plt.savefig( - 'results/plots/fig_ntr_{}_{:.0f}Hz_{}of{}_threshold_{:.2e}V.png'.format( - cfg['trigger_name'], - cfg['target_global_trigger_rate'] / units.Hz, - cfg['number_coincidences'], - cfg['total_number_triggered_channels'], - cfg['final_threshold'])) - plt.show() - plt.close() diff --git a/NuRadioReco/examples/cr_efficiency_analysis/4_air_shower_reconstruction.py b/NuRadioReco/examples/cr_efficiency_analysis/4_air_shower_reconstruction.py deleted file mode 100644 index 09deb30c9..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/4_air_shower_reconstruction.py +++ /dev/null @@ -1,186 +0,0 @@ -import numpy as np -import os, glob -import json, pickle -import helper_cr_eff as hcr -from NuRadioReco.utilities import units -from NuRadioReco.detector.generic_detector import GenericDetector -import NuRadioReco.modules.io.coreas.readCoREASStation -import NuRadioReco.modules.efieldToVoltageConverter -import NuRadioReco.modules.RNO_G.hardwareResponseIncorporator -import NuRadioReco.modules.channelGenericNoiseAdder -import NuRadioReco.modules.channelGalacticNoiseAdder -import NuRadioReco.modules.trigger.envelopeTrigger -import NuRadioReco.modules.trigger.powerIntegration -import NuRadioReco.modules.trigger.highLowThreshold -import NuRadioReco.modules.channelBandPassFilter -import NuRadioReco.modules.eventTypeIdentifier -import NuRadioReco.modules.channelResampler -import NuRadioReco.modules.electricFieldResampler -import NuRadioReco.modules.io.eventWriter -import logging -import argparse - -logger = logging.getLogger() -logger.setLevel(logging.INFO) - -''' -This script reconstructs the air shower (stored in air_shower_sim as hdf5 files) with -the trigger parameters calculated before''' - -parser = argparse.ArgumentParser(description='Run air shower Reconstruction') - -parser.add_argument('--detector_file', type=str, nargs='?', - default='example_data/arianna_station_32.json', - help='choose detector for air shower simulation') -parser.add_argument('--default_station', type=int, nargs='?', default=32, - help='define default station for detector') -parser.add_argument('--triggered_channels', type=list, nargs='?', - default=[0, 1, 2, 3], help='define channels with trigger') -parser.add_argument('--config_file', type=str, nargs='?', - default='config/air_shower/final_config_power_integration_trigger_0Hz_2of3_3.00mV.json', - help='settings from the ntr results') -parser.add_argument('--eventlist', type=str, nargs='?', - default=['example_data/example_data.hdf5'], help='list with event files') -parser.add_argument('--number', type=int, nargs='?', - default=0, help='number of element in eventlist') -parser.add_argument('--output_filename', type=str, nargs='?', - default='output_air_shower_reco/air_shower_reco_', help='begin of output filename') - -args = parser.parse_args() -input_file = args.eventlist[args.number] - -os.makedirs(args.output_filename, exist_ok=True) - -config_file = glob.glob('{}*'.format(args.config_file))[0] -with open(config_file, 'r') as fp: - cfg = json.load(fp) - -logger.info(f"Apply {config_file} on {args.detector_file} with default station {args.default_station} on {input_file}") - -det = GenericDetector(json_filename=args.detector_file, default_station=args.default_station) -det.update(cfg['station_time']) - -station_ids = det.get_station_ids() -station_id = station_ids[0] -channel_ids = det.get_channel_ids(station_id) - -# initialize all modules that are needed for processing -# provide input parameters that are to remain constant during processing -readCoREASStation = NuRadioReco.modules.io.coreas.readCoREASStation.readCoREASStation() -readCoREASStation.begin([input_file], args.default_station, debug=False) -efieldToVoltageConverter = NuRadioReco.modules.efieldToVoltageConverter.efieldToVoltageConverter() -efieldToVoltageConverter.begin(debug=False) -hardwareResponseIncorporator = NuRadioReco.modules.RNO_G.hardwareResponseIncorporator.hardwareResponseIncorporator() -channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() -channelGenericNoiseAdder.begin() -channelGalacticNoiseAdder = NuRadioReco.modules.channelGalacticNoiseAdder.channelGalacticNoiseAdder() -channelGalacticNoiseAdder.begin( - n_side=cfg['galactic_noise_n_side'], - interpolation_frequencies=np.arange(cfg['galactic_noise_interpolation_frequencies_start'], - cfg['galactic_noise_interpolation_frequencies_stop'], - cfg['galactic_noise_interpolation_frequencies_step'])) -if cfg['trigger_name'] == 'high_low': - triggerSimulator = NuRadioReco.modules.trigger.highLowThreshold.triggerSimulator() - triggerSimulator.begin() - -if cfg['trigger_name'] == 'envelope': - triggerSimulator = NuRadioReco.modules.trigger.envelopeTrigger.triggerSimulator() - triggerSimulator.begin() - -if cfg['trigger_name'] == 'power_integration': - triggerSimulator = NuRadioReco.modules.trigger.power_integration.triggerSimulator() - triggerSimulator.begin() - -channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() -channelBandPassFilter.begin() -eventTypeIdentifier = NuRadioReco.modules.eventTypeIdentifier.eventTypeIdentifier() -electricFieldResampler = NuRadioReco.modules.electricFieldResampler.electricFieldResampler() -electricFieldResampler.begin() -channelResampler = NuRadioReco.modules.channelResampler.channelResampler() -channelResampler.begin() -eventWriter = NuRadioReco.modules.io.eventWriter.eventWriter() -eventWriter.begin(args.output_filename - + str(cfg['trigger_name']) + '_' - + str(round(cfg['target_global_trigger_rate'] / units.Hz, 0)) + 'Hz_' - + str(cfg['number_coincidences']) + 'of' - + str(cfg['total_number_triggered_channels']) + '_' - + str(round(cfg['final_threshold']/units.mV, 2)) + 'mV_' - + str(args.number) + '.nur') - -# Loop over all events in file as initialized in readCoRREAS and perform analysis -i = 0 -for evt in readCoREASStation.run(detector=det): - for sta in evt.get_stations(): - #if i == 10: # use this if you want to test something or if you want only 10 position - # break - logger.info("processing event {:d} with id {:d}".format(i, evt.get_id())) - - if cfg['station_time_random']: - sta = hcr.set_random_station_time(sta, cfg['station_time']) - - efieldToVoltageConverter.run(evt, sta, det) - - eventTypeIdentifier.run(evt, sta, "forced", 'cosmic_ray') - channelGenericNoiseAdder.run(evt, sta, det, amplitude=cfg['Vrms_thermal_noise'], - min_freq=cfg['T_noise_min_freq'], max_freq=cfg['T_noise_max_freq'], - type='rayleigh') - channelGalacticNoiseAdder.run(evt, sta, det) - - if cfg['hardware_response']: - hardwareResponseIncorporator.run(evt, sta, det, sim_to_data=True) - - - # in the high_low and power_integration the filter is applied externally - if cfg['trigger_name'] == 'high_low': - channelBandPassFilter.run(evt, sta, det, passband=cfg['passband_trigger'], - filter_type='butter', order=cfg['order_trigger']) - - triggerSimulator.run(evt, sta, det, - threshold_high=cfg['final_threshold'], - threshold_low=-cfg['final_threshold'], - coinc_window=cfg['coinc_window'], - number_concidences=cfg['number_coincidences'], - triggered_channels=args.triggered_channels, - trigger_name='{}_pb_{:.0f}_{:.0f}_tt_{:.2f}'.format( - cfg['trigger_name'], - cfg['passband_trigger'][0] / units.MHz, - cfg['passband_trigger'][1] / units.MHz, - cfg['final_threshold'] / units.mV)) - - if cfg['trigger_name'] == 'envelope': - # The bandpass for the envelope trigger is included in the trigger module, - triggerSimulator.run(evt, sta, det, - passband=cfg['passband_trigger'], - order=cfg['order_trigger'], - number_coincidences=cfg['number_coincidences'], - threshold=cfg['final_threshold'], - coinc_window=cfg['coinc_window'], - triggered_channels=args.triggered_channels, - trigger_name='{}_pb_{:.0f}_{:.0f}_tt_{:.2f}'.format( - cfg['trigger_name'], - cfg['passband_trigger'][0] / units.MHz, - cfg['passband_trigger'][1] / units.MHz, - cfg['final_threshold'] / units.mV)) - - if cfg['trigger_name'] == 'power_integration': - channelBandPassFilter.run(evt, sta, det, passband=cfg['passband_trigger'], - filter_type='butter', order=cfg['order_trigger']) - - triggerSimulator.run(evt, sta, det, - threshold=cfg['final_threshold'], - integration_window=cfg['int_window'], - number_concidences=cfg['number_coincidences'], - triggered_channels=args.triggered_channels, - coinc_window=cfg['coinc_window'], - trigger_name='{}_pb_{:.0f}_{:.0f}_tt_{:.2f}'.format( - cfg['trigger_name'], - cfg['passband_trigger'][0] / units.MHz, - cfg['passband_trigger'][1] / units.MHz, - cfg['final_threshold'] / units.mV)) - - channelResampler.run(evt, sta, det, sampling_rate=cfg['sampling_rate']) - - electricFieldResampler.run(evt, sta, det, sampling_rate=cfg['sampling_rate']) - i += 1 - eventWriter.run(evt, det=det, mode='micro') # here you can change what should be stored in the nur files -nevents = eventWriter.end() diff --git a/NuRadioReco/examples/cr_efficiency_analysis/5_analyze_air_shower_reco.py b/NuRadioReco/examples/cr_efficiency_analysis/5_analyze_air_shower_reco.py deleted file mode 100644 index 83b50dc43..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/5_analyze_air_shower_reco.py +++ /dev/null @@ -1,257 +0,0 @@ -import numpy as np -import os -import glob -import helper_cr_eff as hcr -import json -import NuRadioReco.modules.io.eventReader as eventReader -from NuRadioReco.framework.parameters import stationParameters as stnp -from NuRadioReco.framework.parameters import showerParameters as shp -from NuRadioReco.utilities import units -import argparse -import logging -logger = logging.getLogger() -logger.setLevel(logging.INFO) - -''' -This script sorts the triggered events into different energy, zenith, and distance bins. -Energy means shower energy, -zenith bin means inclination of the shower arrival direction and -distance means distance between shower core and station. -These are important parameters to determine the trigger efficiency. -''' - -parser = argparse.ArgumentParser(description='Nurfile analyser') -parser.add_argument('--config_file', type=str, nargs='?', - default='config/air_shower/final_config_envelope_trigger_1Hz_2of3_22.46mV.json', - help='settings from threshold analysis') -parser.add_argument('--directory', type=str, nargs='?', - default='output_air_shower_reco/', - help='path were results from air shower analysis are stored') -parser.add_argument('--condition', type=str, nargs='?', default='envelope_trigger_0Hz_3of4', - help='string which should be in dict name') -parser.add_argument('--energy_binning', type=list, nargs='?', - default=[16.5, 20, 8], help='energy bins as log() with start, stop, number of bins (np.logspace)') -parser.add_argument('--zenith_binning', type=list, nargs='?', - default=[0, 100, 10], help='zenith bins in deg with start, stop, step (np.arange)') -parser.add_argument('--distance_binning', type=int, nargs='?', - default=[0, 700, 4000], help='distance bin edges as list') -parser.add_argument('--number_of_sta_in_evt', type=int, nargs='?', - default=72, help='number of stations in one event') - -args = parser.parse_args() -energy_binning = args.energy_binning -zenith_binning = args.zenith_binning -distance_binning = args.distance_binning -number_of_sta_in_evt = args.number_of_sta_in_evt - -# the entries of this list are defined in the input argument energy_bins. -# [0] is the start value, [1] is the stop value, [2] is the number of samples generated -energy_bins = np.logspace(energy_binning[0], energy_binning[1], energy_binning[2]) -energy_bins_low = energy_bins[0:-1] -energy_bins_high = energy_bins[1:] -logger.info(f"Use energy bins {energy_bins} eV") - -# the entries of this list are defined in the input argument zenith_bins. -# [0] is the start value, [1] is the stop value, [2] is step size -zenith_bins = np.arange(zenith_binning[0], zenith_binning[1], zenith_binning[2]) * units.deg -zenith_bins_low = zenith_bins[0:-1] -zenith_bins_high = zenith_bins[1:] -logger.info(f"Use zenith bins {zenith_bins/units.deg} deg") - -# the entries of this list are defined in the input argument distance_bins. -distance_bins_low = np.array(distance_binning[0:-1]) -distance_bins_high = np.array(distance_binning[1:]) -logger.info(f"Use distance bins low:{distance_bins_low} to high:{distance_bins_high} m") - -nur_file_list = [] # get non corrupted input files with specified passband -i = 0 -for nur_file in glob.glob('{}*.nur'.format(args.directory)): - if os.path.isfile(nur_file) and str(args.condition) in nur_file: - i += 1 - if os.path.getsize(nur_file) > 0: - nur_file_list.append(nur_file) - -n_files = len(nur_file_list) - -energy = [] -zenith = [] -azimuth = [] -distance = [] -events = [] -trigger_status = [] # trigger status per event station and threshold -trigger_status_weight = [] # trigger status per event station and threshold with weighting -trigger_in_station = [] # name of trigger - -weight = [] -num = 0 - -evtReader = eventReader.eventReader() -evtReader.begin(filename=nur_file, read_detector=True) -det = evtReader.get_detector() -default_station = det.get_station_ids()[0] -for evt in evtReader.run(): # loop over all events, one event is one station - num += 1 - event_id = evt.get_id() - events.append(evt) - det_position = det.get_absolute_position(station_id=default_station) - sta = evt.get_station(station_id=default_station) - sim_station = sta.get_sim_station() - energy.append(sim_station.get_parameter(stnp.cr_energy)) - zenith.append(sim_station.get_parameter(stnp.zenith)) # get zenith for each station - azimuth.append(sim_station.get_parameter(stnp.azimuth)) - current_weight = sim_station.get_simulation_weight() / (units.m**2) - weight.append(current_weight) - trigger_in_station.append(sta.get_triggers()) # get trigger for each station - trigger_status.append(sta.has_triggered()) - trigger_status_weight.append(sta.has_triggered() * current_weight) - - for sho in evt.get_sim_showers(): - core = sho.get_parameter(shp.core) - distance.append(np.sqrt( - ((core[0] - det_position[0])**2) - + (core[1] - det_position[1])**2)) - -trigger_status = np.array(trigger_status) -trigger_status_weight = np.array(trigger_status_weight) -zenith = np.array(zenith) -zenith_deg = zenith / units.deg # this is necessary to avoid mistakes due to decimals -distance = np.array(distance) -energy = np.array(energy) - -# here we reshape the array in a form that the shower parameter are stored once instead one entry for each station. -# Energy and Zenith are shower parameters. -energy_shower = np.array(energy).reshape(int(len(energy)/number_of_sta_in_evt), number_of_sta_in_evt)[:, 0] -zenith_shower = np.array(zenith).reshape(int(len(zenith)/number_of_sta_in_evt), number_of_sta_in_evt)[:, 0] - -# here we calculate the footprint of the shower, e.g. the area which is covered by the shower -footprint_shower = np.sum(np.array(weight).reshape(int(len(weight)/number_of_sta_in_evt), number_of_sta_in_evt), axis=1) - -# here we calculate the area of the footprint where the shower triggers a station -footprint_triggered_area_shower = np.sum(np.array(trigger_status_weight). - reshape(int(len(trigger_status_weight)/number_of_sta_in_evt), - number_of_sta_in_evt), axis=1) - -# here is the trigger status sorted by shower -trigger_status_shower = np.sum(np.array(trigger_status). - reshape(int(len(trigger_status)/number_of_sta_in_evt), - number_of_sta_in_evt), axis=1) - -# here we create empty array which will be filled in the following loop. The first axis contains all parameters in -# the energy bin, the second axis the zenthis bins and the third axis the distance bins - -# number of true triggered trigger in each bin -triggered_trigger_e = np.zeros((len(energy_bins_low), len(zenith_bins_low), len(distance_bins_low))) -# the weight (in this case this is the area) of true triggered for each bin -triggered_trigger_weight_e = np.zeros((len(energy_bins_low), len(zenith_bins_low))) -# events within a bin -masked_events_e = np.zeros((len(energy_bins_low), len(zenith_bins_low), len(distance_bins_low))) -# trigger efficiency in each bin -trigger_efficiency_e = np.zeros((len(energy_bins_low), len(zenith_bins_low), len(distance_bins_low))) -# effective area of each bin (= area inwhich a event within that bin triggers) -trigger_effective_area_e = np.zeros((len(energy_bins_low), len(zenith_bins_low))) -# error of the effective area of each energy bin -trigger_effective_area_err_e = np.zeros((len(energy_bins_low), len(zenith_bins_low))) - - -for dim_0, energy_bin_low, energy_bin_high in zip(range(len(energy_bins_low)), energy_bins_low, energy_bins_high): - mask_e = (energy >= energy_bin_low) & (energy < energy_bin_high) # choose one energy bin - n_events_masked_e = np.sum(mask_e) # number of events in that energy bin - - for dim_1, zenith_bin_low, zenith_bin_high in zip(range(len(zenith_bins_low)), zenith_bins_low, zenith_bins_high): - # choose zenith bin - mask_z = (zenith_deg >= zenith_bin_low/units.deg) & (zenith_deg < zenith_bin_high/units.deg) - # trigger in in one energy bin and one zenith bin (ez) (values depend on the loop) - mask_ez = mask_e & mask_z - masked_trigger_status_ez = trigger_status[mask_ez] - masked_trigger_status_weight_ez = trigger_status_weight[mask_ez] - n_events_masked_ez = np.sum(mask_ez) # number of events in that energy and zenith bin - - # mask ez - # number of triggered true trigger in energy and zentih bin - triggered_trigger_ez = np.sum(masked_trigger_status_ez, axis=0) - triggered_trigger_weight = np.sum(masked_trigger_status_weight_ez, axis=0) - # fraction of events in this energy and zeniths bin that triggered true - trigger_efficiency_ez = triggered_trigger_ez / n_events_masked_ez - - # reshape the array. zenith and energy are the same for all stations in the shower - triggered_trigger_weight_shower = masked_trigger_status_weight_ez.reshape( - int(len(masked_trigger_status_weight_ez) / number_of_sta_in_evt), number_of_sta_in_evt) - - # array with the effective area of all showers in this energy and zenith bin - triggered_trigger_weight_shower_sum = np.sum(triggered_trigger_weight_shower, axis=1) - trigger_effective_area = np.mean(triggered_trigger_weight_shower_sum) - trigger_effective_area_err = np.std(triggered_trigger_weight_shower_sum) - - # set the values of this energy bin in an array for all bins - triggered_trigger_weight_e[dim_0, dim_1] = triggered_trigger_weight - trigger_effective_area_e[dim_0, dim_1] = trigger_effective_area - trigger_effective_area_err_e[dim_0, dim_1] = trigger_effective_area_err - - for dim_2, distance_bin_low, distance_bin_high in zip(range(len(distance_bins_low)), - distance_bins_low, distance_bins_high): - # bin for each event, since distance between shower core and station differs for each station - # choose here, if distances should be in circles or rings, - # so if it should include everything with in or only an interval - mask_d = (distance < distance_bin_high) - mask_ezd = mask_ez & mask_d - masked_trigger_status_ezd = trigger_status[mask_ezd] - masked_trigger_status_weight_ezd = trigger_status_weight[mask_ezd] - n_events_masked_ezd = np.sum(mask_ezd) - - # mask ezd, no a eff - triggered_trigger_ezd = np.sum(masked_trigger_status_ezd, axis=0) - trigger_efficiency_ezd = triggered_trigger_ezd / n_events_masked_ezd - - masked_events_e[dim_0, dim_1, dim_2] = n_events_masked_ezd - triggered_trigger_e[dim_0, dim_1, dim_2] = triggered_trigger_ezd - trigger_efficiency_e[dim_0, dim_1, dim_2] = trigger_efficiency_ezd - -with open(args.config_file, 'r') as fp: - cfg = json.load(fp) - -dic = {} -for key in cfg: - dic[key] = cfg[key] -dic['detector_file'] = [] -dic['default_station'] = default_station -dic['energy_bins_low'] = energy_bins_low -dic['energy_bins_high'] = energy_bins_high -dic['zenith_bins_low'] = zenith_bins_low -dic['zenith_bins_high'] = zenith_bins_high -dic['distance_bins_high'] = distance_bins_high -dic['distance_bins_low'] = distance_bins_low -dic['trigger_effective_area'] = np.nan_to_num(trigger_effective_area_e) -dic['trigger_effective_area_err'] = np.nan_to_num(trigger_effective_area_err_e) -dic['trigger_masked_events'] = np.nan_to_num(masked_events_e) -dic['triggered_trigger'] = np.nan_to_num(triggered_trigger_e) -dic['trigger_efficiency'] = np.nan_to_num(trigger_efficiency_e) -dic['triggered_trigger_weight'] = triggered_trigger_weight_e - -if cfg['hardware_response']: - out_dir = 'results/air_shower_{}_trigger_{:.0f}Hz_{}of{}_{:.2f}mV'.format( - cfg['trigger_name'], - cfg['target_global_trigger_rate'] / units.Hz, - cfg['number_coincidences'], - cfg['total_number_triggered_channels'], - cfg['final_threshold'] / units.millivolt - ) - -else: - out_dir = 'results/air_shower_{}_trigger_{:.0f}Hz_{}of{}_{:.2f}V'.format( - cfg['trigger_name'], - cfg['target_global_trigger_rate'] / units.Hz, - cfg['number_coincidences'], - cfg['total_number_triggered_channels'], - cfg['final_threshold'] - ) - -os.makedirs(out_dir, exist_ok=True) - -output_file = 'dict_air_shower_e{}_z{}_d{}_{}.json'.format( - len(energy_bins_low), - len(zenith_bins_low), - len(distance_bins_low), - max(distance_binning)) - -with open(os.path.join(out_dir, output_file), 'w') as outfile: - json.dump(dic, outfile, cls=hcr.NumpyEncoder, indent=4) diff --git a/NuRadioReco/examples/cr_efficiency_analysis/7_plot_trigger_efficiency.py b/NuRadioReco/examples/cr_efficiency_analysis/7_plot_trigger_efficiency.py deleted file mode 100644 index 7bacf829d..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/7_plot_trigger_efficiency.py +++ /dev/null @@ -1,66 +0,0 @@ -import numpy as np -import json -import matplotlib.pyplot as plt - -filename = 'results/air_shower/dict_air_shower_pb_80_180_e7_z9_d2_4000.json' - -with open(filename, 'r') as fp: - data = json.load(fp) - -distance_bins_low = np.array(data['distance_bins_low']) -energy_bins_low = np.array(data['energy_bins_low']) -zenith_bins_low = np.array(data['zenith_bins_low']) -distance_bins_high = np.array(data['distance_bins_high']) -energy_bins_high = np.array(data['energy_bins_high']) -zenith_bins_high = np.array(data['zenith_bins_high']) - -energy_center = (energy_bins_low + energy_bins_high) / 2 - -trigger_efficiency = np.array(data['trigger_efficiency']) -triggered_trigger = np.array(data['triggered_trigger']) -masked_events = np.array(data['trigger_masked_events']) - - -### choose only some bins -trigger_efficiency = np.nan_to_num(trigger_efficiency[:,0:8,:]) -zenith_bins_low = zenith_bins_low[0:8] -zenith_bins_high = zenith_bins_high[0:8] -triggered_trigger = np.nan_to_num(triggered_trigger[1:3,0:8,:]) -masked_events = np.nan_to_num(masked_events[1:3,0:8,:]) - -triggered_trigger_ez = np.sum(triggered_trigger, axis=2) -masked_events_ez = np.sum(masked_events, axis=2) -trigger_efficiency_ez = triggered_trigger_ez / masked_events_ez - -plt.plot(energy_center[:5], trigger_efficiency[:5, 0, 0], marker='x') -plt.hist(energy_center[:5], bins=np.logspace(15, 19, 9), weights=trigger_efficiency[:5, 0, 0], edgecolor='k', color='lightgrey') -plt.xscale('log') -plt.xlim(1e15, 1e20) -plt.xlabel('Cosmic ray energy [eV]', fontsize=18) -plt.ylabel('Trigger efficiency', fontsize=18) -plt.show() -plt.close() - - -zen = (zenith_bins_low + zenith_bins_high) /2 -energy = energy_bins_low - -int = np.linspace(0, 1, len(distance_bins_low)+1) -for i_energy in range(len(energy)): - plt.plot(zen, trigger_efficiency_ez[i_energy, :], marker='x', label=r'Energy {}, all distances'.format(energy[i_energy]), color='k', linestyle='dashed' ) - for i_dist in range(len(distance_bins_low)): - if i_dist == len(distance_bins_low) - 1: - label = r'd < {} m'.format(distance_bins_high[i_dist]) - plt.plot(zen, trigger_efficiency[i_energy, :, i_dist], marker='x', label=label, color='k') - else: - label = r'd < {} m'.format(distance_bins_high[i_dist]) - plt.plot(zen, trigger_efficiency[i_energy,:,i_dist], marker='x', label=label) - plt.xlabel('Zenith $[^\circ]$', fontsize=18) - plt.ylabel(r'Trigger efficiency', fontsize=18) - plt.title('CR energy {:.0e} $-$ {:.0e}'.format(energy_bins_low[i_energy], energy_bins_high[i_energy])) - plt.xlim(0,90) - plt.ylim(-0.01,1.01) - plt.legend(fontsize=10) - plt.tick_params(axis='both', labelsize=16) - plt.tight_layout() - plt.show() diff --git a/NuRadioReco/examples/cr_efficiency_analysis/LPDA_Southpole.json b/NuRadioReco/examples/cr_efficiency_analysis/LPDA_Southpole.json deleted file mode 100644 index 3d70f3769..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/LPDA_Southpole.json +++ /dev/null @@ -1,30 +0,0 @@ -{ - "_default": {}, - "channels": { - "1": { - "amp_type": "rno_surface", - "ant_comment": "upward facing LPDA", - "ant_orientation_phi": 0.0, - "ant_orientation_theta": 0.0, - "ant_position_x": 0.0, - "ant_position_y": 0.0, - "ant_position_z": -3.0, - "ant_rotation_phi": 180.0, - "ant_rotation_theta": 90.0, - "ant_type": "createLPDA_100MHz", - "channel_id": 1, - "station_id": 101 - } - }, - "positions": {}, - "stations": { - "1": { - "pos_altitude": 0.0, - "pos_easting": 0.0, - "pos_northing": 0, - "pos_site": "southpole", - "station_id": 101, - "station_type": null - } - } -} diff --git a/NuRadioReco/examples/cr_efficiency_analysis/plot_threshold_passband_comparison.py b/NuRadioReco/examples/cr_efficiency_analysis/plot_threshold_passband_comparison.py deleted file mode 100644 index da3d8b6c4..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/plot_threshold_passband_comparison.py +++ /dev/null @@ -1,111 +0,0 @@ -import numpy as np -import json -import os -import matplotlib.pyplot as plt -from NuRadioReco.utilities import units, io_utilities -import argparse - -'''This script plots the results from the analysis for different passbands. It makes only sense if you want to compare -different trigger settings such as passband which are stored in different dict. You can adjust the binning accordingly.''' - -parser = argparse.ArgumentParser(description='plot') -parser.add_argument('n_energy_bins', type=int, nargs='?', default=7, help='number of energy bins in the dict') -parser.add_argument('n_zenith_bins', type=int, nargs='?', default=9, help='number of zenith bins in the dict') -parser.add_argument('n_distance_bins', type=int, nargs='?', default=2, help='number of zenith bins in the dict') -parser.add_argument('n_passbands', type=int, nargs='?', default=1, help='number of different passbands or dict files') - -# choose index of zenith bins for plotting here -zenith_start = int(0 /10) -zenith_stop = int(90 /10) - -args = parser.parse_args() -n_energy_bins = args.n_energy_bins -n_zenith_bins = args.n_zenith_bins -n_distance_bins = args.n_distance_bins -n_passbands = args.n_passbands - -number_coincidences_list = [] -coinc_windows = [] -passband_trigger = np.zeros((n_passbands, 2)) -trigger_threshold = np.zeros((n_passbands)) - -# number of passbands, energy bins, zenith angle bins, distance bins -triggered_trigger = np.zeros((n_passbands, n_energy_bins, n_zenith_bins, n_distance_bins)) -trigger_weight = np.zeros((n_passbands, n_energy_bins, n_zenith_bins, n_distance_bins)) -masked_events = np.zeros((n_passbands, n_energy_bins, n_zenith_bins, n_distance_bins)) -trigger_efficiency = np.zeros((n_passbands, n_energy_bins, n_zenith_bins, n_distance_bins)) -trigger_effective_area = np.zeros((n_passbands, n_energy_bins, n_zenith_bins)) -trigger_effective_area_err = np.zeros((n_passbands, n_energy_bins, n_zenith_bins)) - -pos = 0 -for file in os.listdir('results/air_shower/'): - filename = os.path.join('results/air_shower/', file) - if filename.endswith('.json'): - if os.path.getsize(filename) > 0: - with open(filename, 'r') as fp: - data = json.load(fp) - T_noise = data['T_noise'] - T_noise_min_freq = data['T_noise_min_freq'] - T_noise_max_freq = data['T_noise_max_freq'] - order_trigger = data['order_trigger'] - coinc_window = data['coinc_window'] - number_coincidences = data['number_coincidences'] - - energy_bins_low = data['energy_bins_low'] - energy_bins_high = data['energy_bins_high'] - zenith_bins_low = data['zenith_bins_low'] - zenith_bins_high = data['zenith_bins_high'] - - trigger_threshold[pos] = data['threshold'] - passband_trigger[pos] = data['passband_trigger'] - trigger_efficiency[pos] = data['trigger_efficiency'] - triggered_trigger[pos] = data['triggered_trigger'] - masked_events[pos] = data['trigger_masked_events'] - trigger_effective_area[pos] = data['trigger_effective_area'] - trigger_effective_area_err[pos] = data['trigger_effective_area_err'] - trigger_weight = data['triggered_trigger_weight'] - - pos += 1 - -trigger_efficiency = trigger_efficiency[:, :, zenith_start:zenith_stop] -triggered_trigger = triggered_trigger[:, :, zenith_start:zenith_stop] -masked_events = masked_events[:, :, zenith_start:zenith_stop] -trigger_effective_area = trigger_effective_area[:, :, zenith_start:zenith_stop] -trigger_effective_area_err = trigger_effective_area_err[:, :, zenith_start:zenith_stop] -zenith_bin = np.arange(zenith_start, zenith_stop, 10) - -weight_area = (-(np.cos((zenith_bin+10)*units.deg)**2)/2) - (-(np.cos(zenith_bin*units.deg)**2)/2) #integral of sin*cos is the weight -aeff_zenith = trigger_effective_area * weight_area -total_aeff = np.nansum(aeff_zenith, axis=2) # sum over all zenith bins to get number of cr for each energy bin - -print(aeff_zenith,aeff_zenith.shape, total_aeff) -x_min = np.min(passband_trigger[:, 0])/units.megahertz -x_max = np.max(passband_trigger[:, 0])/units.megahertz -y_min = np.min(passband_trigger[:, 1])/units.megahertz -y_max = np.max(passband_trigger[:, 1])/units.megahertz - -binWidth = 10 -binLength = 10 -x_edges = np.arange(x_min-5, x_max+15, binWidth) -y_edges = np.arange(y_min-5, y_max+15, binLength) - -for count_energy, energy_bin in enumerate(energy_bins_low): - x = passband_trigger[:, 0] / units.megahertz - y = passband_trigger[:, 1] / units.megahertz - z = total_aeff[:,count_energy]/units.km**2 - - counts, xbins, ybins, image = plt.hist2d(x, y, - bins=[x_edges, y_edges], weights=z, - vmin=min(z), vmax=max(z), cmap=plt.cm.jet, cmin = 1e-9) - plt.colorbar(label=r'Effective area [$km^2$]') - plt.title(r'Energy {:.2e}, zenith {}$^\circ$ $-$ {}$^\circ$ '.format(energy_bin, zenith_start*10, zenith_stop*10)) - plt.xlim(x_min - 5, x_max + 5) - plt.ylim(y_min - 5, y_max + 5) - plt.xlabel('Lower cutoff frequency [MHz]') - plt.ylabel('Upper cutoff frequency [MHz]') - plt.legend() - plt.tight_layout() - plt.show() - plt.savefig('results/air_shower/fig_pb_thresh_signal_area_energy_{:.2e}_{}_{}.png'.format( - energy_bin, zenith_start*10, zenith_stop*10)) - plt.close() diff --git a/NuRadioReco/examples/cr_efficiency_analysis/read_me.txt b/NuRadioReco/examples/cr_efficiency_analysis/read_me.txt deleted file mode 100644 index 0212fbacf..000000000 --- a/NuRadioReco/examples/cr_efficiency_analysis/read_me.txt +++ /dev/null @@ -1,32 +0,0 @@ -This directory contains all necessary files to do a cosmic ray analysis for a specific detector, -including trigger settings and air shower reconstruction. - -1. '1_create_config_file.py' -Create a configuration file with all relevant parameters. - -2. '2_calculate_trigger_threshold_for_trigger_rate.py' - or - '2_I_calculate_trigger_rate_for_thresholds.py' and '2_II_merge_output theshold_calulcations.py' -In this step, the threshold for a specific antenna set is calculated. One can choose between a brut force -script which increases the threshold incrementally until the target trigger rate is obtained or a -smarter and faster script which estimates the threshold range and calculates the trigger rate for the given -threshold range. The advantage is, that the hugh number of iterations can be divided into several job and -will be merged with the script 'merge_output theshold_calulcations.py' Even if you only do one job with -'calculate_trigger_rate_for_thresholds.py ' you have to merge the one file afterwards! - -3. '3_plot_trigger_rate_vs_threshold.py' -This script plots the slope of the trigger rate for different thresholds obtained -by 2_calculate_trigger_threshold_for_trigger_rate.py -or '2_I_calculate_trigger_rate_for_thresholds.py' and '2_II_merge_output theshold_calulcations.py' - -here you can also use the extra script plot_threshold_passband_comparison.py - -4. '4_air_shower_reconstruction.py' -This script evaluates the obtained trigger settings on air shower simulations done with corsika and -stores the results in nur files. - -5. '5_anaylze_air_shower_reco.py' -To analyze the performancec of the detector, the results of the air shower reco have to be divided into -different energy, zenith and distance bins. this is done by 'anaylze_air_shower_reco.py' - -6. plot your results with 6_.. and 7_... \ No newline at end of file diff --git a/NuRadioReco/framework/electric_field.py b/NuRadioReco/framework/electric_field.py index e0f3bed48..1e780217e 100644 --- a/NuRadioReco/framework/electric_field.py +++ b/NuRadioReco/framework/electric_field.py @@ -23,7 +23,7 @@ def __init__(self, channel_ids, position=None, Parameters ---------- - channel_ids: array of ints + channel_ids: list of ints the channels ids this electric field is valid for. (For cosmic rays one electric field is typically valid for several channels. For neutrino simulations, we typically diff --git a/NuRadioReco/modules/electricFieldSignalReconstructor.py b/NuRadioReco/modules/electricFieldSignalReconstructor.py index 4d8163f88..305f9321b 100644 --- a/NuRadioReco/modules/electricFieldSignalReconstructor.py +++ b/NuRadioReco/modules/electricFieldSignalReconstructor.py @@ -60,16 +60,16 @@ def run(self, evt, station, det, debug=False): electric_field[efp.signal_time] = signal_time # - low_pos = int(130 * units.ns * electric_field.get_sampling_rate()) - up_pos = int(210 * units.ns * electric_field.get_sampling_rate()) - if(debug): - fig, ax = plt.subplots(1, 1) - sc = ax.scatter(trace_copy[1, low_pos:up_pos], trace_copy[2, low_pos:up_pos], c=electric_field.get_times()[low_pos:up_pos], s=5) - fig.colorbar(sc, ax=ax) - ax.set_aspect('equal') - ax.set_xlabel("eTheta") - ax.set_ylabel("ePhi") - fig.tight_layout() + # low_pos = int(130 * units.ns * electric_field.get_sampling_rate()) + # up_pos = int(210 * units.ns * electric_field.get_sampling_rate()) + # if(debug): + # fig, ax = plt.subplots(1, 1) + # sc = ax.scatter(trace_copy[1, low_pos:up_pos], trace_copy[2, low_pos:up_pos], c=electric_field.get_times()[low_pos:up_pos], s=5) + # fig.colorbar(sc, ax=ax) + # ax.set_aspect('equal') + # ax.set_xlabel("eTheta") + # ax.set_ylabel("ePhi") + # fig.tight_layout() low_pos, up_pos = hp.get_interval(envelope_mag, scale=0.5) v_start = trace_copy[:, signal_time_bin] @@ -99,7 +99,8 @@ def run(self, evt, station, det, debug=False): mask_signal_window = (times > (signal_time - self.__signal_window_pre)) & (times < (signal_time + self.__signal_window_post)) mask_noise_window = np.zeros_like(mask_signal_window, dtype=bool) if(self.__noise_window > 0): - mask_noise_window[int(np.round((-self.__noise_window - 141.) * electric_field.get_sampling_rate())):int(np.round(-141. * electric_field.get_sampling_rate()))] = np.ones(int(np.round(self.__noise_window * electric_field.get_sampling_rate())), dtype=bool) # the last n bins + # set the noise window to the first "self.__noise_window" ns of the trace. If this cuts into the signal window, the noise window is reduced to not overlap with the signal window + mask_noise_window = times < min(times[0] + self.__noise_window, signal_time - self.__signal_window_pre) signal_energy_fluence = trace_utilities.get_electric_field_energy_fluence(trace, times, mask_signal_window, mask_noise_window) dt = times[1] - times[0] @@ -125,13 +126,14 @@ def run(self, evt, station, det, debug=False): electric_field.set_parameter_error(efp.polarization_angle, pol_angle_error) # compute expeted polarization - site = det.get_site(station.get_id()) - exp_efield = hp.get_lorentzforce_vector(electric_field[efp.zenith], electric_field[efp.azimuth], hp.get_magnetic_field_vector(site)) - cs = coordinatesystems.cstrafo(electric_field[efp.zenith], electric_field[efp.azimuth], site=site) - exp_efield_onsky = cs.transform_from_ground_to_onsky(exp_efield) - exp_pol_angle = np.arctan2(np.abs(exp_efield_onsky[2]), np.abs(exp_efield_onsky[1])) - logger.info("expected polarization angle = {:.1f}".format(exp_pol_angle / units.deg)) - electric_field.set_parameter(efp.polarization_angle_expectation, exp_pol_angle) + if det is not None: + site = det.get_site(station.get_id()).lower() + exp_efield = hp.get_lorentzforce_vector(electric_field[efp.zenith], electric_field[efp.azimuth], hp.get_magnetic_field_vector(site)) + cs = coordinatesystems.cstrafo(electric_field[efp.zenith], electric_field[efp.azimuth], site=site) + exp_efield_onsky = cs.transform_from_ground_to_onsky(exp_efield) + exp_pol_angle = np.arctan2(np.abs(exp_efield_onsky[2]), np.abs(exp_efield_onsky[1])) + logger.info("expected polarization angle = {:.1f}".format(exp_pol_angle / units.deg)) + electric_field.set_parameter(efp.polarization_angle_expectation, exp_pol_angle) def end(self): pass diff --git a/NuRadioReco/modules/io/coreas/coreas.py b/NuRadioReco/modules/io/coreas/coreas.py index bb9c0564d..806eafd4c 100644 --- a/NuRadioReco/modules/io/coreas/coreas.py +++ b/NuRadioReco/modules/io/coreas/coreas.py @@ -1,16 +1,22 @@ import numpy as np import matplotlib.pyplot as plt +from matplotlib import cm from radiotools import helper as hp from radiotools import coordinatesystems + from NuRadioReco.utilities import units import NuRadioReco.framework.sim_station +import NuRadioReco.framework.event import NuRadioReco.framework.base_trace import NuRadioReco.framework.electric_field import NuRadioReco.framework.radio_shower -import radiotools.coordinatesystems from NuRadioReco.framework.parameters import stationParameters as stnp from NuRadioReco.framework.parameters import electricFieldParameters as efp from NuRadioReco.framework.parameters import showerParameters as shp +import copy +import h5py +import cr_pulse_interpolator.interpolation_fourier +import cr_pulse_interpolator.signal_interpolation_fourier import logging logger = logging.getLogger('NuRadioReco.coreas') @@ -19,40 +25,511 @@ conversion_fieldstrength_cgs_to_SI = 2.99792458e10 * units.micro * units.volt / units.meter -def get_angles(corsika): +def get_angles(corsika, declination): """ - Converting angles in corsika coordinates to local coordinates + Converting angles in corsika coordinates to local coordinates. + + Corsika positive x-axis points to the magnetic north, NRR coordinates positive x-axis points to the geographic east. + Corsika positive y-axis points to the west, NRR coordinates positive y-axis points to the geographic north. + Corsika z-axis points upwards, NuRadio z-axis points upwards. + + Corsika's zenith angle of a particle trajectory is defined between the particle momentum vector and the negative + z-axis, meaning that the particle is described in the direction where it is going to. The azimuthal angle is + described between the positive x-axis and the horizontal component of the particle momentum vector + (i.e. with respect to the magnetic north) proceeding counterclockwise. + + NRR describes the particle zenith and azimuthal angle in the direction where the particle is coming from. + Therefore the zenith angle is the same, but the azimuthal angle has to be shifted by 180 + 90 degrees. + The north has to be shifted by 90 degrees plus difference between geomagetic and magnetic north. + + + Parameters + ---------- + corsika : hdf5 file object + the open hdf5 file object of the corsika hdf5 file + declination : float + declination of the magnetic field, in internal units + + Returns + ------- + zenith : float + zenith angle + azimuth : float + azimuth angle + magnetic_field_vector : np.ndarray + magnetic field vector """ - zenith = np.deg2rad(corsika['inputs'].attrs["THETAP"][0]) - azimuth = hp.get_normalized_angle(3 * np.pi / 2. + np.deg2rad(corsika['inputs'].attrs["PHIP"][0])) + # TODO: verify sign of declination - Bx, Bz = corsika['inputs'].attrs["MAGNET"] - B_inclination = np.arctan2(Bz, Bx) + zenith = corsika['inputs'].attrs["THETAP"][0] * units.deg + azimuth = hp.get_normalized_angle( + 3 * np.pi / 2. + np.deg2rad(corsika['inputs'].attrs["PHIP"][0]) + declination / units.rad + ) * units.rad - B_strength = (Bx ** 2 + Bz ** 2) ** 0.5 * units.micro * units.tesla + # in CORSIKA convention, the first component points North (y in NRR) and the second component points down (minus z) + By, minBz = corsika['inputs'].attrs["MAGNET"] + B_inclination = np.arctan2(minBz, By) # angle from y-axis towards negative z-axis - # in local coordinates north is + 90 deg - magnetic_field_vector = B_strength * hp.spherical_to_cartesian(np.pi * 0.5 + B_inclination, 0 + np.pi * 0.5) + B_strength = np.sqrt(By ** 2 + minBz ** 2) * units.micro * units.tesla + + # zenith of the magnetic field vector is 90 deg + inclination, as inclination proceeds downwards from horizontal + # azimuth of the magnetic field vector is 90 deg - declination, as declination proceeds clockwise from North + magnetic_field_vector = B_strength * hp.spherical_to_cartesian( + np.pi / 2 + B_inclination, np.pi / 2 - declination / units.rad + ) return zenith, azimuth, magnetic_field_vector +def get_geomagnetic_angle(zenith, azimuth, magnetic_field_vector): + """ + Calculates the angle between the geomagnetic field and the shower axis defined by `zenith` and `azimuth`. + + Parameters + ---------- + zenith : float + zenith angle (in internal units) + azimuth : float + azimuth angle (in internal units) + magnetic_field_vector : np.ndarray + The magnetic field vector in the NRR coordinate system (x points East, y points North, z points up) + + Returns + ------- + geomagnetic_angle : float + geomagnetic angle + """ + shower_axis_vector = hp.spherical_to_cartesian(zenith / units.rad, azimuth / units.rad) + geomagnetic_angle = hp.get_angle(magnetic_field_vector, shower_axis_vector) * units.rad + + return geomagnetic_angle + + +def convert_obs_to_nuradio_efield(observer, zenith, azimuth, magnetic_field_vector): + """ + Converts the electric field from one CoREAS observer to NuRadio units and the on-sky coordinate system. + + The on-sky CS in NRR has basis vectors eR, eTheta, ePhi. + Before going to the on-sky CS, we account for the magnetic field which does not point strictly North. + To get the zenith, azimuth and magnetic field vector, one can use `get_angles()`. + The `observer` array should have the shape (n_samples, 4) with the columns (time, Ey, -Ex, Ez), + where (x, y, z) is the NuRadio CS. + + Parameters + ---------- + observer : np.ndarray + The observer as in the HDF5 file, e.g. list(corsika['CoREAS']['observers'].values())[i]. + zenith : float + zenith angle (in internal units) + azimuth : float + azimuth angle (in internal units) + magnetic_field_vector : np.ndarray + magnetic field vector + + Returns + ------- + efield: np.array (3, n_samples) + Electric field in the on-sky CS (r, theta, phi) + efield_times: np.array (n_samples) + The time values corresponding to the electric field samples + + """ + cs = coordinatesystems.cstrafo( + zenith / units.rad, azimuth / units.rad, + magnetic_field_vector # the magnetic field vector is used to find showerplane, so only direction is important + ) + + efield_times = observer[:, 0] * units.second + efield = np.array([ + observer[:, 2] * -1, # CORSIKA y-axis points West + observer[:, 1], + observer[:, 3] + ]) * conversion_fieldstrength_cgs_to_SI + + # convert coreas efield to NuRadio spherical coordinated eR, eTheta, ePhi (on sky) + efield_geographic = cs.transform_from_magnetic_to_geographic(efield) + efield_on_sky = cs.transform_from_ground_to_onsky(efield_geographic) + + return efield_on_sky, efield_times + + +def convert_obs_positions_to_nuradio_on_ground(observer_pos, zenith, azimuth, magnetic_field_vector): + """ + Convert observer positions from the CORSIKA CS to the NRR ground CS. + + First, the observer position is converted to the NRR coordinate conventions (i.e. x pointing East, + y pointing North, z pointing up). Then, the observer position is corrected for magnetic north + (as CORSIKA only has two components to its magnetic field vector) and put in the geographic CS. + To get the zenith, azimuth and magnetic field vector, one can use the `get_angles()` function. + If multiple observers are to be converted, the `observer` array should have the shape (n_observers, 3). + + Parameters + ---------- + observer_pos : np.ndarray + The observer's position as extracted from the HDF5 file, e.g. corsika['CoREAS']['my_observer'].attrs['position'] + zenith : float + zenith angle (in internal units) + azimuth : float + azimuth angle (in internal units) + magnetic_field_vector : np.ndarray + magnetic field vector + + Returns + ------- + obs_positions_geo: np.ndarray + observer positions in geographic coordinates, shaped as (n_observers, 3). + + """ + cs = coordinatesystems.cstrafo( + zenith / units.rad, azimuth / units.rad, + magnetic_field_vector + ) + + # If single position is given, make sure it has the right shape (3,) -> (1, 3) + if observer_pos.ndim == 1: + observer_pos = observer_pos[np.newaxis, :] + + obs_positions = np.array([ + observer_pos[:, 1] * -1, + observer_pos[:, 0], + observer_pos[:, 2] + ]) * units.cm + + # second to last dimension has to be 3 for the transformation + obs_positions_geo = cs.transform_from_magnetic_to_geographic(obs_positions) + + return obs_positions_geo.T + + +def read_CORSIKA7(input_file, declination=None): + """ + this function reads the corsika hdf5 file and returns a sim_station with all relevent information from the file + + Note that the function assumes the energy has been fixed to a single value. + + Parameters + ---------- + input_file : string + path to the corsika hdf5 file + + Returns + ------- + evt : Event + event object containing the corsika information + """ + if declination is None: + declination = 0 + logger.warning( + "No declination given, assuming 0 degrees. This might need to incorrect electric field polarizations." + ) + + corsika = h5py.File(input_file, "r") + + sampling_rate = 1. / (corsika['CoREAS'].attrs['TimeResolution'] * units.second) + zenith, azimuth, magnetic_field_vector = get_angles(corsika, declination) + energy = corsika['inputs'].attrs["ERANGE"][0] * units.GeV # Assume fixed energy + + sim_shower = make_sim_shower(corsika) + + # The traces are stored in a SimStation + sim_station = NuRadioReco.framework.sim_station.SimStation(0) # set sim station id to 0 + + sim_station.set_parameter(stnp.azimuth, azimuth) + sim_station.set_parameter(stnp.zenith, zenith) + sim_station.set_parameter(stnp.cr_energy, energy) + sim_station.set_magnetic_field_vector(magnetic_field_vector) + + sim_station.set_is_cosmic_ray() + + for j_obs, observer in enumerate(corsika['CoREAS']['observers'].values()): + obs_positions_geo = convert_obs_positions_to_nuradio_on_ground( + observer.attrs['position'], zenith, azimuth, magnetic_field_vector + ) + efield, efield_time = convert_obs_to_nuradio_efield( + observer, zenith, azimuth, magnetic_field_vector + ) + + add_electric_field_to_sim_station(sim_station, [j_obs], efield, efield_time[0], zenith, azimuth, sampling_rate, + efield_position=obs_positions_geo) + + evt = NuRadioReco.framework.event.Event(corsika['inputs'].attrs['RUNNR'], corsika['inputs'].attrs['EVTNR']) + stn = NuRadioReco.framework.station.Station(0) # set station id to 0 + stn.set_sim_station(sim_station) + evt.set_station(stn) + evt.add_sim_shower(sim_shower) + corsika.close() + + return evt + + +def plot_footprint_onsky(sim_station, fig=None, ax=None): + """ + plots the footprint of the observer positions in the (vxB, vxvxB) shower plane + + Parameters + ---------- + sim_station : sim station + simulated station object + """ + obs_positions = [] + zz = [] + for efield in sim_station.get_electric_fields(): + obs_positions.append(efield.get_position_onsky()) + zz.append(np.sum(efield[efp.signal_energy_fluence])) + obs_positions = np.array(obs_positions) + zz = np.array(zz) + if fig is None: + fig, ax = plt.subplots() + ax.set_aspect('equal') + sc = ax.scatter(obs_positions[:, 0], obs_positions[:, 1], c=zz, cmap=cm.gnuplot2_r, marker='o', edgecolors='k') + cbar = fig.colorbar(sc, ax=ax, orientation='vertical') + cbar.set_label('fluence [eV/m^2]') + ax.set_xlabel('vxB [m]') + ax.set_ylabel('vx(vxB) [m]') + return fig, ax + + +def create_sim_station(station_id, evt, weight=None): + """ + Creates an NuRadioReco `SimStation` with the information from an `Event` object created with e.g. read_CORSIKA7(). + + Optionally, station can be assigned a weight. + To add an electric field the function add_electric_field_to_sim_station() has to be used. + + Parameters + ---------- + station_id : int + The id to assign to the new station + evt : Event + event object containing the CoREAS output + weight : float + weight corresponds to area covered by station + + Returns + ------- + sim_station: SimStation + simulated station object + """ + coreas_station = evt.get_station(station_id=0) # read_coreas has only station id 0 + coreas_shower = evt.get_first_sim_shower() + coreas_sim_station = coreas_station.get_sim_station() + + # Make the SimStation and store the parameters extracted from the SimShower + sim_station = NuRadioReco.framework.sim_station.SimStation(station_id) + + sim_station.set_parameter(stnp.azimuth, coreas_shower.get_parameter(shp.azimuth)) + sim_station.set_parameter(stnp.zenith, coreas_shower.get_parameter(shp.zenith)) + + sim_station.set_parameter(stnp.cr_energy, coreas_shower.get_parameter(shp.energy)) + sim_station.set_parameter(stnp.cr_xmax, coreas_shower.get_parameter(shp.shower_maximum)) + + sim_station.set_magnetic_field_vector( + coreas_shower.get_parameter(shp.magnetic_field_vector) + ) + + # Check if the high-level attribute is present in the Event + try: + sim_station.set_parameter(stnp.cr_energy_em, coreas_shower.get_parameter(shp.electromagnetic_energy)) + except KeyError: + global warning_printed_coreas_py + if not warning_printed_coreas_py: + logger.warning( + "No high-level quantities in Event, not setting EM energy, this warning will be only printed once" + ) + warning_printed_coreas_py = True + + if coreas_sim_station.is_cosmic_ray(): + sim_station.set_is_cosmic_ray() + + # Set simulation weight + sim_station.set_simulation_weight(weight) + + return sim_station + + +def make_sim_shower(corsika, declination=0): + """ + Creates an NuRadioReco `SimShower` from a CoREAS HDF5 file. + + Parameters + ---------- + corsika : hdf5 file object + the open hdf5 file object of the corsika hdf5 file + declination : float + + Returns + ------- + sim_shower: SimShower + simulated shower object + """ + zenith, azimuth, magnetic_field_vector = get_angles(corsika, declination) + energy = corsika['inputs'].attrs["ERANGE"][0] * units.GeV # Assume fixed energy + + # Create RadioShower to store simulation parameters in Event + sim_shower = NuRadioReco.framework.radio_shower.RadioShower() + + sim_shower.set_parameter(shp.primary_particle, corsika["inputs"].attrs["PRMPAR"]) + sim_shower.set_parameter(shp.observation_level, corsika["inputs"].attrs["OBSLEV"] * units.cm) + + sim_shower.set_parameter(shp.zenith, zenith) + sim_shower.set_parameter(shp.azimuth, azimuth) + sim_shower.set_parameter(shp.magnetic_field_vector, magnetic_field_vector) + sim_shower.set_parameter(shp.energy, energy) + + sim_shower.set_parameter( + shp.core, np.array([ + corsika['CoREAS'].attrs["CoreCoordinateWest"] * -1, + corsika['CoREAS'].attrs["CoreCoordinateNorth"], + corsika['CoREAS'].attrs["CoreCoordinateVertical"] + ]) * units.cm + ) + sim_shower.set_parameter( + shp.shower_maximum, corsika['CoREAS'].attrs['DepthOfShowerMaximum'] * units.g / units.cm2 + ) + sim_shower.set_parameter( + shp.distance_shower_maximum_geometric, corsika['CoREAS'].attrs["DistanceOfShowerMaximum"] * units.cm + ) + sim_shower.set_parameter( + shp.refractive_index_at_ground, corsika['CoREAS'].attrs["GroundLevelRefractiveIndex"] + ) + sim_shower.set_parameter( + shp.magnetic_field_rotation, corsika['CoREAS'].attrs["RotationAngleForMagfieldDeclination"] * units.degree + ) + + if 'ATMOD' in corsika['inputs'].attrs: # this can be false is left on default or when using GDAS atmosphere + sim_shower.set_parameter(shp.atmospheric_model, corsika["inputs"].attrs["ATMOD"]) + + if 'highlevel' in corsika: + sim_shower.set_parameter(shp.electromagnetic_energy, corsika["highlevel"].attrs["Eem"] * units.eV) + else: + global warning_printed_coreas_py + if not warning_printed_coreas_py: + logger.warning( + "No high-level quantities in HDF5 file, not setting EM energy, this warning will be only printed once") + warning_printed_coreas_py = True + + return sim_shower + + +def create_sim_shower(evt, detector=None, station_id=None): + """ + Create an NuRadioReco `SimShower` from an Event object created with e.g. read_CORSIKA7(), + + the core positions are set such that the detector station is on top of + each coreas observer position + + Parameters + ---------- + evt : Event object + event object containing the CoREAS output, e.g. created with read_CORSIKA7() + detector : detector object + station_id : int + station id of the station relative to which the shower core is given + + Returns + ------- + sim_shower: sim shower + simulated shower object + """ + sim_shower = copy.copy(evt.get_first_sim_shower()) + + efields = evt.get_station(0).get_sim_station().get_electric_fields() + + # We can only set the shower core relative to the station if we know its position + if efields is not None and detector is not None and station_id is not None: + efield_pos = [] + for efield in efields: + efield_pos.append(efield.get_position()) + efield_pos = np.array(efield_pos) + station_position = detector.get_absolute_position(station_id) + core_position = station_position - efield_pos + else: + core_position = np.array([0, 0, 0]) + + core_position[2] = sim_shower.get_parameter(shp.observation_level) # do not alter observation level + sim_shower.set_parameter(shp.core, core_position) + + return sim_shower + + +def add_electric_field_to_sim_station( + sim_station, channel_ids, efield, efield_start_time, zenith, azimuth, sampling_rate, efield_position=None +): + """ + Adds an electric field trace to an existing SimStation, with the provided attributes. + + All variables should be provided in internal units. + + Parameters + ---------- + sim_station : SimStation + The simulated station object, e.g. from make_empty_sim_station() + channel_ids : list of int + Channel IDs to associate to the electric field + efield : np.ndarray + Electric field trace, shaped as (3, n_samples) + efield_start_time : float + Start time of the electric field trace + zenith : float + Zenith angle + azimuth : float + Azimuth angle + sampling_rate : float + Sampling rate of the trace + efield_position : np.ndarray or list of float + Position to associate to the electric field + """ + if type(channel_ids) is not list: + channel_ids = [channel_ids] + + electric_field = NuRadioReco.framework.electric_field.ElectricField(channel_ids, position=efield_position) + + electric_field.set_trace(efield, sampling_rate) + electric_field.set_trace_start_time(efield_start_time) + + electric_field.set_parameter(efp.ray_path_type, 'direct') + electric_field.set_parameter(efp.zenith, zenith) + electric_field.set_parameter(efp.azimuth, azimuth) + + sim_station.add_electric_field(electric_field) + + def calculate_simulation_weights(positions, zenith, azimuth, site='summit', debug=False): - """Calculate weights according to the area that one simulated position represents. + """ + Calculate weights according to the area that one simulated position in readCoreasStation represents. Weights are therefore given in units of area. - Note: The volume of a 2d convex hull is the area.""" + Note: The volume of a 2d convex hull is the area. + + Parameters + ---------- + positions : list + station position with [x, y, z] on ground + zenith : float + zenith angle of the shower + azimuth : float + azimuth angle of the shower + site : str + site of the simulation, default is 'summit' + debug : bool + if true, plots are created for debugging + + Returns + ------- + weights : np.array + weights of the observer position + """ import scipy.spatial as spatial positions = np.array(positions) - cstrafo = radiotools.coordinatesystems.cstrafo(zenith=zenith, azimuth=azimuth, magnetic_field_vector=None, - site=site) - x_trafo_from_shower = cstrafo.transform_from_vxB_vxvxB(station_position=np.array([1, 0, 0])) - y_trafo_from_shower = cstrafo.transform_from_vxB_vxvxB(station_position=np.array([0, 1, 0])) - z_trafo_from_shower = cstrafo.transform_from_vxB_vxvxB(station_position=np.array([0, 0, 1])) + cs = coordinatesystems.cstrafo(zenith=zenith, azimuth=azimuth, magnetic_field_vector=None, + site=site) + x_trafo_from_shower = cs.transform_from_vxB_vxvxB(station_position=np.array([1, 0, 0])) + y_trafo_from_shower = cs.transform_from_vxB_vxvxB(station_position=np.array([0, 1, 0])) + z_trafo_from_shower = cs.transform_from_vxB_vxvxB(station_position=np.array([0, 0, 1])) # voronoi has to be calculated in the shower plane due to symmetry reasons - shower = cstrafo.transform_to_vxB_vxvxB(station_position=positions) + shower = cs.transform_to_vxB_vxvxB(station_position=positions) vor = spatial.Voronoi(shower[:, :2]) # algorithm will find no solution if flat simulation is given in 3d. if debug: @@ -69,19 +546,15 @@ def calculate_simulation_weights(positions, zenith, azimuth, site='summit', debu for p in range(0, weights.shape[0]): # loop over all observer positions vertices_shower_2d = vor.vertices[vor.regions[vor.point_region[p]]] - # x_vertice_ground = x_trafo_from_shower[0] * x_vertice_shower + y_trafo_from_shower[0] * y_vertice_shower + z_trafo_from_shower[0] * z_vertice_shower - # y_vertice_ground = x_trafo_from_shower[1] * x_vertice_shower + y_trafo_from_shower[1] * y_vertice_shower + z_trafo_from_shower[1] * z_vertice_shower - # z_vertice_ground = x_trafo_from_shower[2] * x_vertice_shower + y_trafo_from_shower[2] * y_vertice_shower + z_trafo_from_shower[2] * z_vertice_shower - x_vertice_shower = vertices_shower_2d[:, 0] y_vertice_shower = vertices_shower_2d[:, 1] - z_vertice_shower = -(x_trafo_from_shower[2] * x_vertice_shower + y_trafo_from_shower[2] * y_vertice_shower) / z_trafo_from_shower[2] + z_vertice_shower = -( + x_trafo_from_shower[2] * x_vertice_shower + y_trafo_from_shower[2] * y_vertice_shower + ) / z_trafo_from_shower[2] + + vertices_shower_3d = np.column_stack((x_vertice_shower, y_vertice_shower, z_vertice_shower)) - vertices_shower_3d = [] - for iter in range(len(x_vertice_shower)): - vertices_shower_3d.append([x_vertice_shower[iter], y_vertice_shower[iter], z_vertice_shower[iter]]) - vertices_shower_3d = np.array(vertices_shower_3d) - vertices_ground = cstrafo.transform_from_vxB_vxvxB(station_position=vertices_shower_3d) + vertices_ground = cs.transform_from_vxB_vxvxB(station_position=vertices_shower_3d) n_arms = 8 # mask last observer position of each arm length_shower = np.sqrt(shower[:, 0] ** 2 + shower[:, 1] ** 2) @@ -123,136 +596,387 @@ def calculate_simulation_weights(positions, zenith, azimuth, site='summit', debu return weights -def make_sim_station(station_id, corsika, observer, channel_ids, weight=None): - """ - creates an NuRadioReco sim station from the (interpolated) observer object of the coreas hdf5 file - - Parameters - ---------- - station_id : station id - the id of the station to create - corsika : hdf5 file object - the open hdf5 file object of the corsika hdf5 file - observer : hdf5 observer object - channel_ids : - weight : weight of individual station - weight corresponds to area covered by station - - Returns - ------- - sim_station: sim station - simulated station object +class coreasInterpolator: """ + This class provides an interface to interpolate the electric field traces provided by CoREAS. - zenith, azimuth, magnetic_field_vector = get_angles(corsika) - - if(observer is None): - data = np.zeros((512, 4)) - data[:, 0] = np.arange(0, 512) * units.ns / units.second - else: - data = np.copy(observer) - data[:, 1], data[:, 2] = -observer[:, 2], observer[:, 1] - - # convert to SI units - data[:, 0] *= units.second - data[:, 1] *= conversion_fieldstrength_cgs_to_SI - data[:, 2] *= conversion_fieldstrength_cgs_to_SI - data[:, 3] *= conversion_fieldstrength_cgs_to_SI - - cs = coordinatesystems.cstrafo(zenith, azimuth, magnetic_field_vector=magnetic_field_vector) - efield = cs.transform_from_magnetic_to_geographic(data[:, 1:].T) - efield = cs.transform_from_ground_to_onsky(efield) - - # prepend trace with zeros to not have the pulse directly at the start - n_samples_prepend = efield.shape[1] - efield2 = np.zeros((3, n_samples_prepend + efield.shape[1])) - efield2[0] = np.append(np.zeros(n_samples_prepend), efield[0]) - efield2[1] = np.append(np.zeros(n_samples_prepend), efield[1]) - efield2[2] = np.append(np.zeros(n_samples_prepend), efield[2]) - - sampling_rate = 1. / (corsika['CoREAS'].attrs['TimeResolution'] * units.second) - sim_station = NuRadioReco.framework.sim_station.SimStation(station_id) - electric_field = NuRadioReco.framework.electric_field.ElectricField(channel_ids) - electric_field.set_trace(efield2, sampling_rate) - electric_field.set_trace_start_time(data[0, 0]) - electric_field.set_parameter(efp.ray_path_type, 'direct') - electric_field.set_parameter(efp.zenith, zenith) - electric_field.set_parameter(efp.azimuth, azimuth) - sim_station.add_electric_field(electric_field) - sim_station.set_parameter(stnp.azimuth, azimuth) - sim_station.set_parameter(stnp.zenith, zenith) - energy = corsika['inputs'].attrs["ERANGE"][0] * units.GeV - sim_station.set_parameter(stnp.cr_energy, energy) - sim_station.set_magnetic_field_vector(magnetic_field_vector) - sim_station.set_parameter(stnp.cr_xmax, corsika['CoREAS'].attrs['DepthOfShowerMaximum']) - try: - sim_station.set_parameter(stnp.cr_energy_em, corsika["highlevel"].attrs["Eem"]) - except: - global warning_printed_coreas_py - if(not warning_printed_coreas_py): - logger.warning("No high-level quantities in HDF5 file, not setting EM energy, this warning will be only printed once") - warning_printed_coreas_py = True - sim_station.set_is_cosmic_ray() - sim_station.set_simulation_weight(weight) - return sim_station - - -def make_sim_shower(corsika, observer=None, detector=None, station_id=None): - """ - creates an NuRadioReco sim shower from the coreas hdf5 file + The interpolation method is based on the Fourier interpolation method described in + Corstanje et al. (2023), JINST 18 P09005. The implementation lives in a separate package + called cr-pulse-inerpolator (this package is a optional dependency of NuRadioReco). Parameters ---------- - corsika : hdf5 file object - the open hdf5 file object of the corsika hdf5 file - observer : hdf5 observer object - detector : detector object - station_id : station id of the station relativ to which the shower core is given - - Returns - ------- - sim_shower: sim shower - simulated shower object + corsika_evt : Event + An Event object containing the CoREAS output, from read_CORSIKA7() """ - sim_shower = NuRadioReco.framework.radio_shower.RadioShower() - zenith, azimuth, magnetic_field_vector = get_angles(corsika) - sim_shower.set_parameter(shp.zenith, zenith) - sim_shower.set_parameter(shp.azimuth, azimuth) - sim_shower.set_parameter(shp.magnetic_field_vector, magnetic_field_vector) - - energy = corsika['inputs'].attrs["ERANGE"][0] * units.GeV - sim_shower.set_parameter(shp.energy, energy) - # We can only set the shower core relative to the station if we know its position - if observer is not None and detector is not None and station_id is not None: - station_position = detector.get_absolute_position(station_id) - position = observer.attrs['position'] - cs = coordinatesystems.cstrafo(zenith, azimuth, magnetic_field_vector=magnetic_field_vector) - observer_position = np.zeros(3) - observer_position[0], observer_position[1], observer_position[2] = -position[1] * units.cm, position[0] * units.cm, position[2] * units.cm - observer_position = cs.transform_from_magnetic_to_geographic(observer_position) - core_position = (-observer_position + station_position) - core_position[2] = 0 - sim_shower.set_parameter(shp.core, core_position) - - sim_shower.set_parameter(shp.shower_maximum, corsika['CoREAS'].attrs['DepthOfShowerMaximum'] * units.g / units.cm2) - sim_shower.set_parameter(shp.refractive_index_at_ground, corsika['CoREAS'].attrs["GroundLevelRefractiveIndex"]) - sim_shower.set_parameter(shp.magnetic_field_rotation, - corsika['CoREAS'].attrs["RotationAngleForMagfieldDeclination"] * units.degree) - sim_shower.set_parameter(shp.distance_shower_maximum_geometric, - corsika['CoREAS'].attrs["DistanceOfShowerMaximum"] * units.cm) - - sim_shower.set_parameter(shp.observation_level, corsika["inputs"].attrs["OBSLEV"] * units.cm) - sim_shower.set_parameter(shp.primary_particle, corsika["inputs"].attrs["PRMPAR"]) - if 'ATMOD' in corsika['inputs'].attrs.keys(): - sim_shower.set_parameter(shp.atmospheric_model, corsika["inputs"].attrs["ATMOD"]) - - try: - sim_shower.set_parameter(shp.electromagnetic_energy, corsika["highlevel"].attrs["Eem"] * units.eV) - except: - global warning_printed_coreas_py - if(not warning_printed_coreas_py): - logger.warning("No high-level quantities in HDF5 file, not setting EM energy, this warning will be only printed once") - warning_printed_coreas_py = True - - return sim_shower + def __init__(self, corsika_evt): + # These are set by self.initialize_star_shape() + self.sampling_rate = None + self.electric_field_on_sky = None + self.efield_times = None + + self.obs_positions_ground = None + self.obs_positions_showerplane = None + + # Store the SimStation and SimShower objects + self.sim_station = corsika_evt.get_station(0).get_sim_station() + self.shower = corsika_evt.get_first_sim_shower() # there should only be one simulated shower + self.cs = self.shower.get_coordinatesystem() + + # Flags to check whether interpolator is initialized + self.star_shape_initialized = False + self.efield_interpolator_initialized = False + self.fluence_interpolator_initialized = False + + # Interpolator objects + self.interp_lowfreq = None + self.interp_highfreq = None + self.efield_interpolator = None + self.fluence_interpolator = None + + self.initialize_star_shape() + + logger.info( + f'Initialised star shape pattern for interpolation. ' + f'The shower arrives from zenith={self.zenith / units.deg:.1f}deg, ' + f'azimuth={self.azimuth / units.deg:.1f}deg ' + f'The starshape has radius {self.starshape_radius:.0f}m in the shower plane ' + f'and {self.starshape_radius_ground:.0f}m on ground. ' + ) + + def initialize_star_shape(self): + """ + Initializes the star shape pattern for interpolation, e.g. creates the arrays with the observer positions + in the shower plane and the electric field. + """ + obs_positions = [] + electric_field_on_sky = [] + efield_times = [] + + for j_obs, efield in enumerate(self.sim_station.get_electric_fields()): + obs_positions.append(efield.get_position()) + electric_field_on_sky.append(efield.get_trace().T) + efield_times.append(efield.get_times()) + + # shape: (n_observers, n_samples, (eR, eTheta, ePhi)) + self.electric_field_on_sky = np.array(electric_field_on_sky) + self.efield_times = np.array(efield_times) + self.sampling_rate = 1. / (self.efield_times[0][1] - self.efield_times[0][0]) + + self.obs_positions_ground = np.array(obs_positions) # (n_observers, 3) + self.obs_positions_showerplane = self.cs.transform_to_vxB_vxvxB( + self.obs_positions_ground, core=self.shower.get_parameter(shp.core) + ) + + self.star_shape_initialized = True + + @property + def zenith(self): + return self.shower.get_parameter(shp.zenith) + + @property + def azimuth(self): + return self.shower.get_parameter(shp.azimuth) + + @property + def magnetic_field_vector(self): + return self.shower.get_parameter(shp.magnetic_field_vector) + + @property + def starshape_radius(self): + """ + returns the maximal radius of the star shape pattern in the shower plane + """ + if not self.star_shape_initialized: + logger.error('The interpolator was not initialized, call initialize_star_shape first') + return None + else: + return np.max(np.linalg.norm(self.obs_positions_showerplane[:, :-1], axis=-1)) + + @property + def starshape_radius_ground(self): + """ + returns the maximal radius of the star shape pattern on ground + """ + if not self.star_shape_initialized: + logger.error('The interpolator was not initialized, call initialize_star_shape first') + return None + else: + return np.max(np.linalg.norm(self.obs_positions_ground[:, :-1], axis=-1)) + + def get_empty_efield(self): + """ + Get an array of zeros in the shape of the electric field on the sky + """ + if not self.star_shape_initialized: + logger.error('The interpolator was not initialized, call initialize_star_shape first') + return None + else: + return np.zeros_like(self.electric_field_on_sky[0, :, :]) + + def initialize_efield_interpolator(self, interp_lowfreq, interp_highfreq): + """ + Initialise the efield interpolator. + + The efield will be interpolated in the shower plane for + geometrical reasons. If the geomagnetic angle is smaller than 15deg, no interpolator object is returned. + + Parameters + ---------- + interp_lowfreq : float + Lower frequency for the bandpass filter in interpolation (in internal units) + interp_highfreq : float + Upper frequency for the bandpass filter in interpolation (in internal units) + + Returns + ------- + efield_interpolator : interpolator object + + """ + self.interp_lowfreq = interp_lowfreq + self.interp_highfreq = interp_highfreq + + geomagnetic_angle = get_geomagnetic_angle(self.zenith, self.azimuth, self.magnetic_field_vector) + + if geomagnetic_angle < 15 * units.deg: + logger.warning( + f'The geomagnetic angle is {geomagnetic_angle / units.deg:.2f} deg, ' + f'which is smaller than 15deg, which is the lower limit for the signal interpolation. ' + f'The closest observer is used instead.') + self.efield_interpolator = -1 + else: + logger.info( + f'electric field interpolation with lowfreq {interp_lowfreq / units.MHz} MHz ' + f'and highfreq {interp_highfreq / units.MHz} MHz') + + self.efield_interpolator = cr_pulse_interpolator.signal_interpolation_fourier.interp2d_signal( + self.obs_positions_showerplane[:, 0], + self.obs_positions_showerplane[:, 1], + self.electric_field_on_sky, # TODO: this should be a rotated version for showers from zenith or N-S + signals_start_times=self.efield_times[:, 0] / units.s, + lowfreq=(interp_lowfreq - 0.01) / units.MHz, + highfreq=(interp_highfreq + 0.01) / units.MHz, + sampling_period=1 / self.sampling_rate / units.s, # interpolator wants sampling period in seconds + phase_method="phasor", + radial_method='cubic', + upsample_factor=5, + coherency_cutoff_threshold=0.9, + ignore_cutoff_freq_in_timing=False, + verbose=False + ) + + self.efield_interpolator_initialized = True + + return self.efield_interpolator + + def initialize_fluence_interpolator(self, quantity=efp.signal_energy_fluence): + """ + Initialise fluence interpolator. + + Parameters + ---------- + quantity : electric field parameter + quantity to interpolate, e.g. efp.signal_energy_fluence + The quantity needs to be available as parameter in the electric field object!!! You might + need to run the electricFieldSignalReconstructor. The default is efp.signal_energy_fluence. + + Returns + ------- + fluence_interpolator : interpolator object + """ + fluence_per_position = [ + np.sum(efield[quantity]) for efield in self.sim_station.get_electric_fields() + ] # the fluence is calculated per polarization, so we need to sum them up + + logger.info(f'fluence interpolation') + self.fluence_interpolator = cr_pulse_interpolator.interpolation_fourier.interp2d_fourier( + self.obs_positions_showerplane[:, 0], + self.obs_positions_showerplane[:, 1], + fluence_per_position + ) + self.fluence_interpolator_initialized = True + + return self.fluence_interpolator + + def plot_footprint_fluence(self, radius=300): + """ + plots the interpolated values of the fluence in the shower plane + + Parameters + ---------- + radius : float + radius around shower core which should be plotted + + Returns + ------- + fig : figure object + + ax : axis object + """ + + # Make color plot of f(x, y), using a meshgrid + ti = np.linspace(-radius, radius, 500) + XI, YI = np.meshgrid(ti, ti) + + # Get interpolated values at each grid point, calling the instance of interp2d_fourier + ZI = self.fluence_interpolator(XI, YI) + + # And plot it + maxp = np.max(ZI) + fig, ax = plt.subplots(figsize=(8, 6)) + ax.pcolor(XI, YI, ZI, vmax=maxp, vmin=0, cmap=cm.gnuplot2_r) + mm = cm.ScalarMappable(cmap=cm.gnuplot2_r) + mm.set_array([0.0, maxp]) + cbar = plt.colorbar(mm, ax=ax) + cbar.set_label(r'energy fluence [eV/m^2]', fontsize=14) + ax.set_xlabel(r'$\vec{v} \times \vec{B} [m]', fontsize=16) + ax.set_ylabel(r'$\vec{v} \times (\vec{v} \times \vec{B})$ [m]', fontsize=16) + ax.set_xlim(-radius, radius) + ax.set_ylim(-radius, radius) + ax.set_aspect('equal') + return fig, ax + + def get_position_showerplane(self, position_on_ground): + """ + Transform the position of the antenna on ground to the shower plane. + + Parameters + ---------- + position_on_ground : np.ndarray + Position of the antenna on ground + This value can either be a 2D array (x, y) or a 3D array (x, y, z). If the z-coordinate is missing, the + z-coordinate is automatically set to the observation level of the simulation. + """ + core = self.shower.get_parameter(shp.core) + + if len(position_on_ground) == 2: + position_on_ground = np.append(position_on_ground, core[2]) + logger.info( + f"The antenna position is given in 2D, assuming the antenna is on the ground. " + f"The z-coordinate is set to the observation level {core[2] / units.m:.2f}m" + ) + elif position_on_ground[2] != core[2]: + logger.warning( + "The antenna position is not in the same z-plane as the core position. " + "This behaviour is not tested, so only proceed if you know what you are doing." + ) + + antenna_pos_showerplane = self.cs.transform_to_vxB_vxvxB(position_on_ground, core=core) + + return antenna_pos_showerplane + + def get_interp_efield_value(self, position_on_ground): + """ + Calculate the interpolated electric field given an antenna position on ground. + + For the interpolation, the pulse will be projected in the shower plane. + If the geomagnetic angle is smaller than 15deg, the electric field of the closest observer position is returned. + + Parameters + ---------- + position_on_ground : np.ndarray + Position of the antenna on ground + This value can either be a 2D array (x, y) or a 3D array (x, y, z). If the z-coordinate is missing, the + z-coordinate is automatically set to the observation level of the simulation. + + Returns + ------- + efield_interp : float + Interpolated efield value + trace_start_time : float + Start time of the trace + """ + logger.debug( + f"Getting interpolated efield for antenna position {position_on_ground} on ground" + ) + + antenna_pos_showerplane = self.get_position_showerplane(position_on_ground) + logger.debug(f"The antenna position in shower plane is {antenna_pos_showerplane}") + + # calculate distance between core position at(0,0) and antenna positions in shower plane + dcore_showerplane = np.linalg.norm(antenna_pos_showerplane[:-1]) + + # interpolate electric field at antenna position in shower plane which are inside star pattern + if dcore_showerplane > self.starshape_radius: + efield_interp = self.get_empty_efield() + trace_start_time = None + logger.debug( + f'Antenna position with distance {dcore_showerplane:.2f} to core is outside of starshape ' + f'with radius {self.starshape_radius:.2f}. Setting efield to zero and trace_start_time to None' + ) + elif self.efield_interpolator == -1: + efield_interp = self.get_closest_observer_efield(antenna_pos_showerplane) + trace_start_time = None + else: + efield_interp, trace_start_time, _, _ = self.efield_interpolator( + antenna_pos_showerplane[0], antenna_pos_showerplane[1], + lowfreq=self.interp_lowfreq / units.MHz, + highfreq=self.interp_highfreq / units.MHz, + filter_up_to_cutoff=False, + account_for_timing=True, + pulse_centered=True, + full_output=True) + + return efield_interp, trace_start_time + + def get_interp_fluence_value(self, position_on_ground): + """ + Calculate the interpolated fluence for a given position on the ground. + + Parameters + ---------- + position_on_ground : np.ndarray + Position of the antenna on ground + This value can either be a 2D array (x, y) or a 3D array (x, y, z). If the z-coordinate is missing, the + z-coordinate is automatically set to the observation level of the simulation. + + Returns + ------- + fluence_interp : float + interpolated fluence value + """ + logger.debug( + f"Getting interpolated efield for antenna position {position_on_ground} on ground" + ) + + antenna_pos_showerplane = self.get_position_showerplane(position_on_ground) + logger.debug(f"The antenna position in shower plane is {antenna_pos_showerplane}") + + # calculate distance between core position at(0,0) and antenna positions in shower plane + dcore_showerplane = np.linalg.norm(antenna_pos_showerplane[:-1]) + + # interpolate electric field at antenna position in shower plane which are inside star pattern + if dcore_showerplane > self.starshape_radius: + fluence_interp = 0 + logger.debug( + f'Antenna position with distance {dcore_showerplane:.2f} to core is outside of star pattern ' + f'with radius {self.starshape_radius:.2f}, setting fluence to 0') + else: + fluence_interp = self.fluence_interpolator(antenna_pos_showerplane[0], antenna_pos_showerplane[1]) + + return fluence_interp + + def get_closest_observer_efield(self, antenna_pos_showerplane): + """ + Returns the electric field of the closest observer position for an antenna position in the shower plane. + + Parameters + ---------- + antenna_pos_showerplane : np.ndarray + antenna position in the shower plane + + Returns + ------- + efield : array of floats + electric field + + """ + distances = np.linalg.norm(antenna_pos_showerplane[:2] - self.obs_positions_showerplane[:, :2], axis=1) + index = np.argmin(distances) + efield = self.electric_field_on_sky[index, :, :] + logger.debug( + f'Returning the electric field of the closest observer position, ' + f'which is {distances[index] / units.m:.2f}m away from the antenna' + ) + return efield diff --git a/NuRadioReco/modules/io/coreas/readCoREAS.py b/NuRadioReco/modules/io/coreas/readCoREAS.py deleted file mode 100644 index 61ee8ece3..000000000 --- a/NuRadioReco/modules/io/coreas/readCoREAS.py +++ /dev/null @@ -1,205 +0,0 @@ -from NuRadioReco.modules.base.module import register_run -import h5py -import NuRadioReco.framework.event -import NuRadioReco.framework.station -import NuRadioReco.framework.radio_shower -from radiotools import coordinatesystems as cstrafo -from NuRadioReco.modules.io.coreas import coreas -from NuRadioReco.utilities import units -import numpy as np -import numpy.random -import logging -import time -import os - - -class readCoREAS: - - def __init__(self): - self.__t = 0 - self.__t_event_structure = 0 - self.__t_per_event = 0 - self.__input_files = None - self.__station_id = None - self.__n_cores = None - self.__max_distace = None - self.__current_input_file = None - self.__random_generator = None - self.logger = logging.getLogger('NuRadioReco.readCoREAS') - - def begin(self, input_files, station_id, n_cores=10, max_distance=2 * units.km, seed=None): - """ - begin method - - initialize readCoREAS module - - Parameters - ---------- - input_files: input files - list of coreas hdf5 files - station_id: station id - id number of the station - n_cores: number of cores (integer) - the number of random core positions to generate for each input file - max_distance: radius of random cores (double or None) - if None: max distance is set to the maximum ground distance of the - star pattern simulation - seed: int (default: None) - Seed for the random number generation. If None is passed, no seed is set - """ - self.__input_files = input_files - self.__station_id = station_id - self.__n_cores = n_cores - self.__max_distace = max_distance - self.__current_input_file = 0 - - self.__random_generator = numpy.random.RandomState(seed) - - @register_run() - def run(self, detector, output_mode=0): - """ - Read in a random sample of stations from a CoREAS file. - A number of random positions is selected within a certain radius. - For each position the closest observer is selected and a simulated - event is created for that observer. - - Parameters - ---------- - detector: Detector object - Detector description of the detector that shall be simulated - output_mode: integer (default 0) - - * 0: only the event object is returned - * 1: the function reuturns the event object, the current inputfilename, the distance between the choosen station and the requested core position, - and the area in which the core positions are randomly distributed - - - """ - while (self.__current_input_file < len(self.__input_files)): - t = time.time() - t_per_event = time.time() - filesize = os.path.getsize(self.__input_files[self.__current_input_file]) - if(filesize < 18456 * 2): # based on the observation that a file with such a small filesize is corrupt - self.logger.warning("file {} seems to be corrupt, skipping to next file".format(self.__input_files[self.__current_input_file])) - self.__current_input_file += 1 - continue - corsika = h5py.File(self.__input_files[self.__current_input_file], "r") - self.logger.info( - "using coreas simulation {} with E={:2g} theta = {:.0f}".format( - self.__input_files[self.__current_input_file], - corsika['inputs'].attrs["ERANGE"][0] * units.GeV, - corsika['inputs'].attrs["THETAP"][0] - ) - ) - positions = [] - for i, observer in enumerate(corsika['CoREAS']['observers'].values()): - position = observer.attrs['position'] - positions.append(np.array([-position[1], position[0], 0]) * units.cm) - self.logger.debug("({:.0f}, {:.0f})".format(position[0], position[1])) - positions = np.array(positions) - - max_distance = self.__max_distace - if(max_distance is None): - max_distance = np.max(np.abs(positions[:, 0:2])) - area = np.pi * max_distance ** 2 - - if(output_mode == 0): - n_cores = self.__n_cores * 100 # for output mode 1 we want always n_cores in star pattern. Therefore we generate more core positions to be able to select n_cores in the star pattern afterwards - elif(output_mode == 1): - n_cores = self.__n_cores - else: - raise ValueError('output mode {} not defined.'.format(output_mode)) - theta = self.__random_generator.rand(n_cores) * 2 * np.pi - r = (self.__random_generator.rand(n_cores)) ** 0.5 * max_distance - cores = np.array([r * np.cos(theta), r * np.sin(theta), np.zeros(n_cores)]).T - - zenith, azimuth, magnetic_field_vector = coreas.get_angles(corsika) - cs = cstrafo.cstrafo(zenith, azimuth, magnetic_field_vector) - positions_vBvvB = cs.transform_from_magnetic_to_geographic(positions.T) - positions_vBvvB = cs.transform_to_vxB_vxvxB(positions_vBvvB).T - dd = (positions_vBvvB[:, 0] ** 2 + positions_vBvvB[:, 1] ** 2) ** 0.5 - ddmax = dd.max() - self.logger.info("star shape from: {} - {}".format(-dd.max(), dd.max())) - - cores_vBvvB = cs.transform_from_magnetic_to_geographic(cores.T) - cores_vBvvB = cs.transform_to_vxB_vxvxB(cores_vBvvB).T - dcores = (cores_vBvvB[:, 0] ** 2 + cores_vBvvB[:, 1] ** 2) ** 0.5 - mask_cores_in_starpattern = dcores <= ddmax - - if((not np.sum(mask_cores_in_starpattern)) and (output_mode == 1)): # handle special case of no core position being generated within star pattern - observer = corsika['CoREAS']['observers'].values()[0] - - evt = NuRadioReco.framework.event.Event(corsika['inputs'].attrs['RUNNR'], corsika['inputs'].attrs['EVTNR']) # create empty event - station = NuRadioReco.framework.station.Station(self.__station_id) - sim_station = coreas.make_sim_station(self.__station_id, corsika, observer, detector.get_channel_ids(self.__station_id)) - - station.set_sim_station(sim_station) - evt.set_station(station) - yield evt, self.__current_input_file, None, area - - cores_to_iterate = cores_vBvvB[mask_cores_in_starpattern] - if(output_mode == 0): # select first n_cores that are in star pattern - if(np.sum(mask_cores_in_starpattern) < self.__n_cores): - self.logger.warning("only {0} cores contained in star pattern, returning {0} cores instead of {1} cores that were requested".format(np.sum(mask_cores_in_starpattern), self.__n_cores)) - else: - cores_to_iterate = cores_vBvvB[mask_cores_in_starpattern][:self.__n_cores] - - self.__t_per_event += time.time() - t_per_event - self.__t += time.time() - t - - for iCore, core in enumerate(cores_to_iterate): - t = time.time() - # check if out of bounds - - distances = np.linalg.norm(core[:2] - positions_vBvvB[:, :2], axis=1) - index = np.argmin(distances) - distance = distances[index] - key = list(corsika['CoREAS']['observers'].keys())[index] - self.logger.info( - "generating core at ground ({:.0f}, {:.0f}), vBvvB({:.0f}, {:.0f}), nearest simulated station is {:.0f}m away at ground ({:.0f}, {:.0f}), vBvvB({:.0f}, {:.0f})".format( - cores[iCore][0], - cores[iCore][1], - core[0], - core[1], - distance / units.m, - positions[index][0], - positions[index][1], - positions_vBvvB[index][0], - positions_vBvvB[index][1] - ) - ) - t_event_structure = time.time() - observer = corsika['CoREAS']['observers'].get(key) - - evt = NuRadioReco.framework.event.Event(self.__current_input_file, iCore) # create empty event - station = NuRadioReco.framework.station.Station(self.__station_id) - channel_ids = detector.get_channel_ids(self.__station_id) - sim_station = coreas.make_sim_station(self.__station_id, corsika, observer, channel_ids) - station.set_sim_station(sim_station) - evt.set_station(station) - sim_shower = coreas.make_sim_shower(corsika, observer, detector, self.__station_id) - evt.add_sim_shower(sim_shower) - rd_shower = NuRadioReco.framework.radio_shower.RadioShower(station_ids=[station.get_id()]) - evt.add_shower(rd_shower) - if(output_mode == 0): - self.__t += time.time() - t - self.__t_event_structure += time.time() - t_event_structure - yield evt - elif(output_mode == 1): - self.__t += time.time() - t - self.__t_event_structure += time.time() - t_event_structure - yield evt, self.__current_input_file, distance, area - else: - self.logger.debug("output mode > 1 not implemented") - raise NotImplementedError - - self.__current_input_file += 1 - - def end(self): - from datetime import timedelta - self.logger.setLevel(logging.INFO) - dt = timedelta(seconds=self.__t) - self.logger.info("total time used by this module is {}".format(dt)) - self.logger.info("\tcreate event structure {}".format(timedelta(seconds=self.__t_event_structure))) - self.logger.info("per event {}".format(timedelta(seconds=self.__t_per_event))) - return dt diff --git a/NuRadioReco/modules/io/coreas/readCoREASDetector.py b/NuRadioReco/modules/io/coreas/readCoREASDetector.py new file mode 100644 index 000000000..e9bde9704 --- /dev/null +++ b/NuRadioReco/modules/io/coreas/readCoREASDetector.py @@ -0,0 +1,245 @@ +from datetime import timedelta +import logging +import os +import time +import numpy as np +from NuRadioReco.modules.base.module import register_run +import NuRadioReco.framework.event +import NuRadioReco.framework.station +import NuRadioReco.framework.radio_shower +from NuRadioReco.framework.parameters import showerParameters as shp +from NuRadioReco.framework.parameters import electricFieldParameters as efp +import NuRadioReco.modules.io.coreas +from NuRadioReco.modules.io.coreas import coreas +from NuRadioReco.utilities import units +from NuRadioReco.utilities.signal_processing import half_hann_window +from collections import defaultdict + +conversion_fieldstrength_cgs_to_SI = 2.99792458e10 * units.micro * units.volt / units.meter + + +def get_random_core_positions(xmin, xmax, ymin, ymax, n_cores, seed=None): + """ + Generate random core positions within a rectangle + + Parameters + ---------- + xmin: float + minimum x value + xmax: float + maximum x value + ymin: float + minimum y value + ymax: float + maximum y value + n_cores: int + number of cores to generate + seed: int + seed for the random number generator + + Returns + ------- + cores: array (n_cores, 3) + array containing the core positions + """ + random_generator = np.random.RandomState(seed) + + # generate core positions randomly within a rectangle + cores = np.array([random_generator.uniform(xmin, xmax, n_cores), + random_generator.uniform(ymin, ymax, n_cores), + np.zeros(n_cores)]).T + return cores + + +def apply_hanning(efields): + """ + Apply a half-Hann window to the electric field in the time domain. This smoothens the edges + to avoid ringing effects. + + Parameters + ---------- + efield in time domain: array (n_samples, n_polarizations) + + Returns + ------- + smoothed_efield: array (n_samples, n_polarizations) + """ + + if efields is None: + return None + else: + smoothed_trace = np.zeros_like(efields) + filter = half_hann_window(efields.shape[0], half_percent=0.1) + for pol in range(efields.shape[1]): + smoothed_trace[:, pol] = efields[:, pol] * filter + return smoothed_trace + + +def select_channels_per_station(det, station_id, requested_channel_ids): + """ + Returns a defaultdict object containing the requested channel ids that are in the given station. + This dict contains the channel group ids as keys with lists of channel ids as values. + + Parameters + ---------- + det : DetectorBase + The detector object that contains the station + station_id : int + The station id to select channels from + requested_channel_ids : list of int + List of requested channel ids + + Returns + ------- + channel_ids : defaultdict + Dictionary with channel group ids as keys and lists of channel ids as values + """ + channel_ids = defaultdict(list) + for channel_id in requested_channel_ids: + if channel_id in det.get_channel_ids(station_id): + channel_group_id = det.get_channel_group_id(station_id, channel_id) + channel_ids[channel_group_id].append(channel_id) + return channel_ids + + +class readCoREASDetector: + """ + Use this as default when reading CoREAS files and combining them with a detector. + + This module reads the electric fields of a CoREAS file with a star shaped pattern of observers. + The electric field is then interpolated at the positions of the antennas or stations of a detector. + If the angle between magnetic field and shower direction are below about 15 deg, + the interpolation is no longer reliable and the closest observer is used instead. + """ + + def __init__(self): + self.__t = 0 + self.__t_event_structure = 0 + self.__t_per_event = 0 + self.__corsika_evt = None + self.__interp_lowfreq = None + self.__interp_highfreq = None + + self.coreas_interpolator = None + + self.logger = logging.getLogger('NuRadioReco.readCoREASDetector') + + def begin(self, input_file, interp_lowfreq=30 * units.MHz, interp_highfreq=1000 * units.MHz, + log_level=logging.INFO): + """ + begin method, initialize readCoREASDetector module + + Parameters + ---------- + input_file: str + coreas hdf5 file + interp_lowfreq: float, default=30 * units.MHz + lower frequency for the bandpass filter in interpolation, + should be broader than the sensitivity band of the detector + interp_highfreq: float, default=1000 * units.MHz + higher frequency for the bandpass filter in interpolation, + should be broader than the sensitivity band of the detector + log_level: default=logging.INFO + log level for the logger + """ + self.logger.setLevel(log_level) + + filesize = os.path.getsize(input_file) + if filesize < 18456 * 2: # based on the observation that a file with such a small filesize is corrupt + self.logger.warning("file {} seems to be corrupt".format(input_file)) + + self.__corsika_evt = coreas.read_CORSIKA7(input_file) + self.logger.info( + f"using coreas simulation {input_file} with " + f"E={self.__corsika_evt.get_first_sim_shower().get_parameter(shp.energy):.2g}eV, " + f"zenith angle = {self.__corsika_evt.get_first_sim_shower().get_parameter(shp.zenith) / units.deg:.0f}deg" + ) + + self.__interp_lowfreq = interp_lowfreq + self.__interp_highfreq = interp_highfreq + + self.coreas_interpolator = coreas.coreasInterpolator(self.__corsika_evt) + self.coreas_interpolator.initialize_efield_interpolator(self.__interp_lowfreq, self.__interp_highfreq) + + @register_run() + def run(self, detector, core_position_list, selected_station_ids=None, selected_channel_ids=None): + """ + run method, get interpolated electric fields for the given detector and core positions and set them in the event. + The trace is smoothed with a half-Hann window to avoid ringing effects. When using short traces, this might have + a significant effect on the result. + + Parameters + ---------- + detector: Detector object + Detector description of the detector that shall be simulated + core_position_list: list of array (n_cores, 3) + list of core positions in the format [[x1, y1, z1], [x2, y2, z2], ...] + selected_station_ids: list of int + list of station ids to simulate + selected_channel_ids: list of int + list of channel ids to simulate + """ + + if selected_station_ids is None: + selected_station_ids = detector.get_station_ids() + logging.info(f"using all station ids in detector description: {selected_station_ids}") + else: + logging.info(f"using selected station ids: {selected_station_ids}") + + t = time.time() + t_per_event = time.time() + self.__t_per_event += time.time() - t_per_event + self.__t += time.time() - t + + for iCore, core in enumerate(core_position_list): + t = time.time() + evt = NuRadioReco.framework.event.Event(self.__corsika_evt.get_run_number(), iCore) # create empty event + sim_shower = self.__corsika_evt.get_first_sim_shower() + sim_shower.set_parameter(shp.core, core) + evt.add_sim_shower(sim_shower) + # rd_shower = NuRadioReco.framework.radio_shower.RadioShower(station_ids=selected_station_ids) + # evt.add_shower(rd_shower) + corsika_sim_stn = self.__corsika_evt.get_station(0).get_sim_station() + for station_id in selected_station_ids: + station = NuRadioReco.framework.station.Station(station_id) + sim_station = NuRadioReco.framework.sim_station.SimStation(station_id) + for key, value in corsika_sim_stn.get_parameters().items(): + sim_station.set_parameter(key, value) # copy relevant sim_station parameters over + sim_station.set_magnetic_field_vector(corsika_sim_stn.get_magnetic_field_vector()) + sim_station.set_is_cosmic_ray() + + det_station_position = detector.get_absolute_position(station_id) + channel_ids_in_station = detector.get_channel_ids(station_id) + if selected_channel_ids is None: + selected_channel_ids = channel_ids_in_station + channel_ids_dict = select_channels_per_station(detector, station_id, selected_channel_ids) + for ch_g_ids in channel_ids_dict.keys(): + channel_ids_for_group_id = channel_ids_dict[ch_g_ids] + antenna_position_rel = detector.get_relative_position(station_id, channel_ids_for_group_id[0]) + antenna_position = det_station_position + antenna_position_rel + res_efield, res_trace_start_time = self.coreas_interpolator.get_interp_efield_value(antenna_position - core) + smooth_res_efield = apply_hanning(res_efield) + if smooth_res_efield is None: + smooth_res_efield = self.coreas_interpolator.get_empty_efield() + electric_field = NuRadioReco.framework.electric_field.ElectricField(channel_ids_for_group_id) + electric_field.set_trace(smooth_res_efield.T, self.coreas_interpolator.get_sampling_rate()) + electric_field.set_trace_start_time(res_trace_start_time) + electric_field.set_parameter(efp.ray_path_type, 'direct') + electric_field.set_parameter(efp.zenith, sim_shower[shp.zenith]) + electric_field.set_parameter(efp.azimuth, sim_shower[shp.azimuth]) + sim_station.add_electric_field(electric_field) + station.set_sim_station(sim_station) + # distance_to_core = np.linalg.norm(det_station_position[:-1] - core[:-1]) + # station.set_parameter(stnp.distance_to_core, distance_to_core) + evt.set_station(station) + + self.__t += time.time() - t + yield evt + + def end(self): + self.logger.setLevel(logging.INFO) + dt = timedelta(seconds=self.__t) + self.logger.info("total time used by this module is {}".format(dt)) + self.logger.info("\tcreate event structure {}".format(timedelta(seconds=self.__t_event_structure))) + self.logger.info("per event {}".format(timedelta(seconds=self.__t_per_event))) + return dt diff --git a/NuRadioReco/modules/io/coreas/readCoREASShower.py b/NuRadioReco/modules/io/coreas/readCoREASShower.py index f72347bbf..58e01c2cf 100644 --- a/NuRadioReco/modules/io/coreas/readCoREASShower.py +++ b/NuRadioReco/modules/io/coreas/readCoREASShower.py @@ -2,12 +2,8 @@ import NuRadioReco.framework.station from NuRadioReco.framework.parameters import showerParameters as shp from NuRadioReco.modules.io.coreas import coreas -from NuRadioReco.utilities import units -from radiotools import coordinatesystems -import h5py import numpy as np import time -import re import os import logging logger = logging.getLogger('NuRadioReco.coreas.readCoREASShower') @@ -78,39 +74,35 @@ def run(self): t = time.time() t_per_event = time.time() - filesize = os.path.getsize( - self.__input_files[self.__current_input_file]) + filename = self.__input_files[self.__current_input_file] + filesize = os.path.getsize(filename) if filesize < 18456 * 2: # based on the observation that a file with such a small filesize is corrupt logger.warning( - "file {} seems to be corrupt, skipping to next file".format( - self.__input_files[self.__current_input_file] - ) - ) + "file {} seems to be corrupt, skipping to next file".format(filename)) self.__current_input_file += 1 continue logger.info('Reading %s ...' % self.__input_files[self.__current_input_file]) - corsika = h5py.File( - self.__input_files[self.__current_input_file], "r") - logger.info("using coreas simulation {} with E={:2g} theta = {:.0f}".format( - self.__input_files[self.__current_input_file], corsika['inputs'].attrs["ERANGE"][0] * units.GeV, - corsika['inputs'].attrs["THETAP"][0])) - - f_coreas = corsika["CoREAS"] + corsika_evt = coreas.read_CORSIKA7(self.__input_files[self.__current_input_file]) if self.__ascending_run_and_event_number: evt = NuRadioReco.framework.event.Event(self.__ascending_run_and_event_number, self.__ascending_run_and_event_number) self.__ascending_run_and_event_number += 1 else: - evt = NuRadioReco.framework.event.Event( - corsika['inputs'].attrs['RUNNR'], corsika['inputs'].attrs['EVTNR']) - - evt.__event_time = f_coreas.attrs["GPSSecs"] + evt = NuRadioReco.framework.event.Event(corsika_evt.get_run_number, corsika_evt.get_id) # create sim shower, no core is set since no external detector description is given +<<<<<<< HEAD + sim_shower = coreas.create_sim_shower(corsika_evt) + evt.__event_time = sim_shower.get_parameter(shp.coreas_GPSSecs) + + sim_shower.set_parameter(shp.core, np.array([0, 0, sim_shower.get_parameter(shp.core_coordinate_vertical)])) # set core to 0,0 and vertical coordinate + evt.add_sim_shower(sim_shower) + +======= sim_shower = coreas.make_sim_shower(corsika) sim_shower.set_parameter(shp.core, np.array( [-f_coreas.attrs["CoreCoordinateWest"], f_coreas.attrs["CoreCoordinateNorth"], f_coreas.attrs["CoreCoordinateVertical"]] @@ -122,12 +114,44 @@ def run(self): sim_shower.get_parameter(shp.zenith), sim_shower.get_parameter(shp.azimuth), magnetic_field_vector=sim_shower.get_parameter(shp.magnetic_field_vector)) +>>>>>>> develop # add simulated pulses as sim station - for idx, (name, observer) in enumerate(f_coreas['observers'].items()): - # returns proper station id if possible - station_id = antenna_id(name, idx) + corsika_efields = corsika_evt.get_station(0).get_sim_station().get_electric_fields() + for station_id, corsika_efield in enumerate(corsika_efields): station = NuRadioReco.framework.station.Station(station_id) +<<<<<<< HEAD + sim_station = coreas.create_sim_station(station_id, corsika_evt) + efield_trace = corsika_efield.get_trace() + efield_sampling_rate = efield_trace.get_sampling_rate() + efield_times = corsika_efield.get_times() + + prepend_zeros = True # prepend zeros to not have the pulse directly at the start, heritage from old code + if prepend_zeros: + n_samples_prepend = efield_trace.shape[1] + efield_cor = np.zeros((3, n_samples_prepend + efield_trace.shape[1])) + efield_cor[0] = np.append(np.zeros(n_samples_prepend), efield_trace[0]) + efield_cor[1] = np.append(np.zeros(n_samples_prepend), efield_trace[1]) + efield_cor[2] = np.append(np.zeros(n_samples_prepend), efield_trace[2]) + + efield_times_cor = np.arange(0, n_samples_prepend + efield_trace.shape[1]) / efield_sampling_rate + + else: + efield_cor = efield_trace + efield_times_cor = efield_times + + if self.__det is None: + channel_ids = [0] + else: + if self.__det.has_station(station_id): + channel_ids = self.__det.get_channel_ids(station_id) + else: + channel_ids = self.__det.get_channel_ids(self.__det.get_reference_station_ids()[0]) + + coreas.add_electric_field_to_sim_station(sim_station, channel_ids, efield_cor, efield_times_cor[0], + sim_shower.get_parameter(shp.zenith), + sim_shower.get_parameter(shp.azimuth), efield_sampling_rate) +======= if self.__det is None: sim_station = coreas.make_sim_station( station_id, corsika, observer, channel_ids=[0, 1]) @@ -136,27 +160,32 @@ def run(self): station_id, corsika, observer, channel_ids=self.__det.get_channel_ids(self.__det.get_default_station_id())) +>>>>>>> develop station.set_sim_station(sim_station) evt.set_station(station) if self.__det is not None: +<<<<<<< HEAD + efield_pos = corsika_efield.get_position() +======= position = observer.attrs['position'] antenna_position = np.array([-position[1], position[0], position[2]]) * units.cm antenna_position = cs.transform_from_magnetic_to_geographic(antenna_position) +>>>>>>> develop if not self.__det.has_station(station_id): self.__det.add_generic_station({ 'station_id': station_id, - 'pos_easting': antenna_position[0], - 'pos_northing': antenna_position[1], - 'pos_altitude': antenna_position[2], + 'pos_easting': efield_pos[0], + 'pos_northing': efield_pos[1], + 'pos_altitude': efield_pos[2], 'reference_station': self.__det.get_reference_station_ids()[0] }) else: self.__det.add_station_properties_for_event({ - 'pos_easting': antenna_position[0], - 'pos_northing': antenna_position[1], - 'pos_altitude': antenna_position[2] + 'pos_easting': efield_pos[0], + 'pos_northing': efield_pos[1], + 'pos_altitude': efield_pos[2] }, station_id, evt.get_run_number(), evt.get_id()) self.__t_per_event += time.time() - t_per_event @@ -177,17 +206,4 @@ def end(self): timedelta(seconds=self.__t_event_structure))) logger.info("per event {}".format( timedelta(seconds=self.__t_per_event))) - return dt - - -def antenna_id(antenna_name, default_id): - """ - This function parses the antenna names given in a CoREAS simulation and tries to find an ID - It can be extended to other name patterns - """ - - if re.match("AERA_", antenna_name): - new_id = int(antenna_name.strip("AERA_")) - return new_id - else: - return default_id + return dt \ No newline at end of file diff --git a/NuRadioReco/modules/io/coreas/readCoREASStation.py b/NuRadioReco/modules/io/coreas/readCoREASStation.py index f0bdc8e16..5db4d0b43 100644 --- a/NuRadioReco/modules/io/coreas/readCoREASStation.py +++ b/NuRadioReco/modules/io/coreas/readCoREASStation.py @@ -1,15 +1,20 @@ from NuRadioReco.modules.base.module import register_run -import h5py import numpy as np import NuRadioReco.framework.event import NuRadioReco.framework.station from NuRadioReco.modules.io.coreas import coreas -from NuRadioReco.utilities import units +from NuRadioReco.framework.parameters import stationParameters as stnp import logging logger = logging.getLogger('NuRadioReco.coreas.readCoREASStation') class readCoREASStation: + """ + Reads in CoREAS simulations and creates simulated events for each observer, i.e., a new event for each simulated observer. + + This module is useful for studies of individual electric fields, e.g., to detemine how well the energy fluence + can be reconstructed as a function of singal-to-noise ratio, or to study the polarization reconstruction. + """ def begin(self, input_files, station_id, debug=False): """ @@ -44,49 +49,62 @@ def run(self, detector): """ for input_file in self.__input_files: self.__current_event = 0 - with h5py.File(input_file, "r") as corsika: - if "highlevel" not in corsika.keys() or list(corsika["highlevel"].values()) == []: - logger.warning(" No highlevel quantities in simulated hdf5 files, weights will be taken from station position") - positions = [] - for observer in corsika['CoREAS']['observers'].values(): - position = observer.attrs['position'] - positions.append(np.array([-position[1], position[0], 0]) * units.cm) - positions = np.array(positions) - zenith, azimuth, magnetic_field_vector = coreas.get_angles(corsika) - weights = coreas.calculate_simulation_weights(positions, zenith, azimuth, debug=self.__debug) - - if self.__debug: - import matplotlib.pyplot as plt - fig, ax = plt.subplots() - im = ax.scatter(positions[:, 0], positions[:, 1], c=weights) - fig.colorbar(im, ax=ax).set_label(label=r'Area $[m^2]$') - plt.xlabel('East [m]') - plt.ylabel('West [m]') - plt.title('Final weighting') - plt.gca().set_aspect('equal') - plt.show() + + corsika_evt = coreas.read_CORSIKA7(input_file) + coreas_sim_station = corsika_evt.get_station(0).get_sim_station() + corsika_efields = coreas_sim_station.get_electric_fields() + + efield_pos = [] + for corsika_efield in corsika_efields: + efield_pos.append(corsika_efield.get_position()) + efield_pos = np.array(efield_pos) + + weights = coreas.calculate_simulation_weights(efield_pos, coreas_sim_station.get_parameter(stnp.zenith), coreas_sim_station.get_parameter(stnp.azimuth), debug=self.__debug) + if self.__debug: + import matplotlib.pyplot as plt + fig, ax = plt.subplots() + im = ax.scatter(efield_pos[:, 0], efield_pos[:, 1], c=weights) + fig.colorbar(im, ax=ax).set_label(label=r'Area $[m^2]$') + plt.xlabel('East [m]') + plt.ylabel('West [m]') + plt.title('Final weighting') + plt.gca().set_aspect('equal') + plt.show() + + # make one event for each observer with differen core position (set in create_sim_shower) + for i, corsika_efield in enumerate(corsika_efields): + evt = NuRadioReco.framework.event.Event(self.__current_input_file, self.__current_event) # create empty event + station = NuRadioReco.framework.station.Station(self.__station_id) + sim_station = coreas.create_sim_station(self.__station_id, corsika_evt, weights[i]) + + channel_ids = detector.get_channel_ids(self.__station_id) + efield_trace = corsika_efield.get_trace() + efield_sampling_rate = corsika_efield.get_sampling_rate() + efield_times = corsika_efield.get_times() + + prepend_zeros = True # prepend zeros to not have the pulse directly at the start, heritage from old code + if prepend_zeros: + n_samples_prepend = efield_trace.shape[1] + efield_cor = np.zeros((3, n_samples_prepend + efield_trace.shape[1])) + efield_cor[0] = np.append(np.zeros(n_samples_prepend), efield_trace[0]) + efield_cor[1] = np.append(np.zeros(n_samples_prepend), efield_trace[1]) + efield_cor[2] = np.append(np.zeros(n_samples_prepend), efield_trace[2]) + + efield_times_cor = np.arange(0, n_samples_prepend + efield_trace.shape[1]) / efield_sampling_rate else: - positions = list(corsika["highlevel"].values())[0]["antenna_position"] - zenith, azimuth, magnetic_field_vector = coreas.get_angles(corsika) - weights = coreas.calculate_simulation_weights(positions, zenith, azimuth) - - for i, (name, observer) in enumerate(corsika['CoREAS']['observers'].items()): - evt = NuRadioReco.framework.event.Event(self.__current_input_file, self.__current_event) # create empty event - station = NuRadioReco.framework.station.Station(self.__station_id) - sim_station = coreas.make_sim_station( - self.__station_id, - corsika, - observer, - detector.get_channel_ids(self.__station_id), - weights[i] - ) - station.set_sim_station(sim_station) - sim_shower = coreas.make_sim_shower(corsika, observer, detector, self.__station_id) - evt.add_sim_shower(sim_shower) - evt.set_station(station) - self.__current_event += 1 - yield evt + efield_cor = efield_trace + efield_times_cor = efield_times + + coreas.add_electric_field_to_sim_station(sim_station, channel_ids, efield_cor, efield_times_cor[0], + sim_station.get_parameter(stnp.zenith), + sim_station.get_parameter(stnp.azimuth), efield_sampling_rate) + station.set_sim_station(sim_station) + sim_shower = coreas.create_sim_shower(corsika_evt, detector, self.__station_id) + evt.add_sim_shower(sim_shower) + evt.set_station(station) + self.__current_event += 1 + yield evt self.__current_input_file += 1 def end(self): diff --git a/NuRadioReco/modules/io/coreas/readCoREASStationGrid.py b/NuRadioReco/modules/io/coreas/readCoREASStationGrid.py deleted file mode 100644 index 142b20fe7..000000000 --- a/NuRadioReco/modules/io/coreas/readCoREASStationGrid.py +++ /dev/null @@ -1,211 +0,0 @@ -from NuRadioReco.modules.base.module import register_run -import h5py -import NuRadioReco.framework.event -import NuRadioReco.framework.station -import NuRadioReco.framework.radio_shower -from radiotools import coordinatesystems as cstrafo -from NuRadioReco.framework.parameters import showerParameters as shp -from NuRadioReco.modules.io.coreas import coreas -from NuRadioReco.utilities import units -import numpy as np -import numpy.random -import logging -import time -import os - - -class readCoREAS: - """ - coreas input module for fixed grid of stations. - This module distributes core positions randomly within a user defined area and calculates the electric field - at the detector positions as specified in the detector description by choosing the closest antenna - of the star shape pattern simulation - """ - - def __init__(self): - self.__t = 0 - self.__t_event_structure = 0 - self.__t_per_event = 0 - self.__input_files = None - self.__station_id = None - self.__n_cores = None - self.__max_distace = None - self.__current_input_file = None - self.__random_generator = None - self.logger = logging.getLogger('NuRadioReco.readCoREAS') - - def begin(self, input_files, xmin, xmax, ymin, ymax, n_cores=10, seed=None, log_level=logging.NOTSET): - """ - begin method - - initialize readCoREAS module - - Parameters - ---------- - input_files: input files - list of coreas hdf5 files - xmin: float - minimum x coordinate of the area in which core positions are distributed - xmax: float - maximum x coordinate of the area in which core positions are distributed - ymin: float - minimum y coordinate of the area in which core positions are distributed - ynax: float - maximum y coordinate of the area in which core positions are distributed - n_cores: number of cores (integer) - the number of random core positions to generate for each input file - seed: int (default: None) - Seed for the random number generation. If None is passed, no seed is set - """ - self.__input_files = input_files - self.__n_cores = n_cores - self.__current_input_file = 0 - self.__area = [xmin, xmax, ymin, ymax] - - self.__random_generator = numpy.random.RandomState(seed) - self.logger.setLevel(log_level) - - @register_run() - def run(self, detector, output_mode=0): - """ - Read in a random sample of stations from a CoREAS file. - For each position the closest observer is selected and a simulated - event is created for that observer. - - Parameters - ---------- - detector: Detector object - Detector description of the detector that shall be simulated - output_mode: integer (default 0) - - * 0: only the event object is returned - * 1: the function reuturns the event object, the current inputfilename, the distance between the choosen station and the requested core position, - and the area in which the core positions are randomly distributed - - - """ - while (self.__current_input_file < len(self.__input_files)): - t = time.time() - t_per_event = time.time() - filesize = os.path.getsize(self.__input_files[self.__current_input_file]) - if(filesize < 18456 * 2): # based on the observation that a file with such a small filesize is corrupt - self.logger.warning("file {} seems to be corrupt, skipping to next file".format(self.__input_files[self.__current_input_file])) - self.__current_input_file += 1 - continue - try: - corsika = h5py.File(self.__input_files[self.__current_input_file], "r") - except OSError as error: - self.logger.error(f"error opening file number {self.__current_input_file}: {self.__input_files[self.__current_input_file]}") - self.logger.error(error) - self.__current_input_file += 1 - continue - - self.logger.info( - "using coreas simulation {} with E={:2g} theta = {:.0f}".format( - self.__input_files[self.__current_input_file], - corsika['inputs'].attrs["ERANGE"][0] * units.GeV, - corsika['inputs'].attrs["THETAP"][0] - ) - ) - positions = [] - for i, observer in enumerate(corsika['CoREAS']['observers'].values()): - position = observer.attrs['position'] - positions.append(np.array([-position[1], position[0], 0]) * units.cm) -# self.logger.debug("({:.0f}, {:.0f})".format(positions[i][0], positions[i][1])) - positions = np.array(positions) - - zenith, azimuth, magnetic_field_vector = coreas.get_angles(corsika) - cs = cstrafo.cstrafo(zenith, azimuth, magnetic_field_vector) - positions_vBvvB = cs.transform_from_magnetic_to_geographic(positions.T) - positions_vBvvB = cs.transform_to_vxB_vxvxB(positions_vBvvB).T -# for i, pos in enumerate(positions_vBvvB): -# self.logger.debug("star shape") -# self.logger.debug("({:.0f}, {:.0f}); ({:.0f}, {:.0f})".format(positions[i, 0], positions[i, 1], pos[0], pos[1])) - - dd = (positions_vBvvB[:, 0] ** 2 + positions_vBvvB[:, 1] ** 2) ** 0.5 - ddmax = dd.max() - self.logger.info("star shape from: {} - {}".format(-dd.max(), dd.max())) - - # generate core positions randomly within a rectangle - cores = np.array([self.__random_generator.uniform(self.__area[0], self.__area[1], self.__n_cores), - self.__random_generator.uniform(self.__area[2], self.__area[3], self.__n_cores), - np.zeros(self.__n_cores)]).T - - self.__t_per_event += time.time() - t_per_event - self.__t += time.time() - t - - station_ids = detector.get_station_ids() - for iCore, core in enumerate(cores): - t = time.time() - evt = NuRadioReco.framework.event.Event(self.__current_input_file, iCore) # create empty event - sim_shower = coreas.make_sim_shower(corsika) - sim_shower.set_parameter(shp.core, core) - evt.add_sim_shower(sim_shower) - rd_shower = NuRadioReco.framework.radio_shower.RadioShower(station_ids=station_ids) - evt.add_shower(rd_shower) - for station_id in station_ids: - # convert into vxvxB frame to calculate closests simulated station to detecor station - det_station_position = detector.get_absolute_position(station_id) - det_station_position[2] = 0 - core_rel_to_station = core - det_station_position - # core_rel_to_station_vBvvB = cs.transform_from_magnetic_to_geographic(core_rel_to_station) - core_rel_to_station_vBvvB = cs.transform_to_vxB_vxvxB(core_rel_to_station) - dcore = (core_rel_to_station_vBvvB[0] ** 2 + core_rel_to_station_vBvvB[1] ** 2) ** 0.5 -# print(f"{core_rel_to_station}, {core_rel_to_station_vBvvB} -> {dcore}") - if(dcore > ddmax): - # station is outside of the star shape pattern, create empty station - station = NuRadioReco.framework.station.Station(station_id) - channel_ids = detector.get_channel_ids(station_id) - sim_station = coreas.make_sim_station(station_id, corsika, None, channel_ids) - station.set_sim_station(sim_station) - evt.set_station(station) - self.logger.debug(f"station {station_id} is outside of star shape, channel_ids {channel_ids}") - else: - distances = np.linalg.norm(core_rel_to_station_vBvvB[:2] - positions_vBvvB[:,:2], axis=1) - index = np.argmin(distances) - distance = distances[index] - key = list(corsika['CoREAS']['observers'].keys())[index] - self.logger.debug( - "generating core at ground ({:.0f}, {:.0f}), rel to station ({:.0f}, {:.0f}) vBvvB({:.0f}, {:.0f}), nearest simulated station is {:.0f}m away at ground ({:.0f}, {:.0f}), vBvvB({:.0f}, {:.0f})".format( - cores[iCore][0], - cores[iCore][1], - core_rel_to_station[0], - core_rel_to_station[1], - core_rel_to_station_vBvvB[0], - core_rel_to_station_vBvvB[1], - distance / units.m, - positions[index][0], - positions[index][1], - positions_vBvvB[index][0], - positions_vBvvB[index][1] - ) - ) - t_event_structure = time.time() - observer = corsika['CoREAS']['observers'].get(key) - - station = NuRadioReco.framework.station.Station(station_id) - channel_ids = detector.get_channel_ids(station_id) - sim_station = coreas.make_sim_station(station_id, corsika, observer, channel_ids) - station.set_sim_station(sim_station) - evt.set_station(station) - if(output_mode == 0): - self.__t += time.time() - t - yield evt - elif(output_mode == 1): - self.__t += time.time() - t - self.__t_event_structure += time.time() - t_event_structure - yield evt, self.__current_input_file - else: - self.logger.debug("output mode > 1 not implemented") - raise NotImplementedError - - self.__current_input_file += 1 - - def end(self): - from datetime import timedelta - self.logger.setLevel(logging.INFO) - dt = timedelta(seconds=self.__t) - self.logger.info("total time used by this module is {}".format(dt)) - self.logger.info("\tcreate event structure {}".format(timedelta(seconds=self.__t_event_structure))) - self.logger.info("per event {}".format(timedelta(seconds=self.__t_per_event))) - return dt diff --git a/NuRadioReco/modules/voltageToEfieldConverterPerChannelGroup.py b/NuRadioReco/modules/voltageToEfieldConverterPerChannelGroup.py new file mode 100644 index 000000000..c5995a80b --- /dev/null +++ b/NuRadioReco/modules/voltageToEfieldConverterPerChannelGroup.py @@ -0,0 +1,107 @@ +from NuRadioReco.modules.base.module import register_run +import numpy as np +from NuRadioReco.utilities import trace_utilities +from NuRadioReco.detector import antennapattern +from NuRadioReco.framework.parameters import stationParameters as stnp +from NuRadioReco.framework.parameters import electricFieldParameters as efp +from NuRadioReco.framework import electric_field as ef +from NuRadioReco.modules.io.coreas.readCoREASDetector import select_channels_per_station +from NuRadioReco.modules.voltageToEfieldConverter import get_array_of_channels +import logging + + +class voltageToEfieldConverterPerChannelGroup: + """ + This module is intended to reconstruct the electric field from dual-polarized antennas, + i.e., two antennas with orthogonal polarizations combined in one mechanical structure. + This is the typical case for air-shower detectors such as Auger and LOFAR. + + Converts voltage trace to electric field per channel group. + """ + + def __init__(self): + self.logger = logging.getLogger('NuRadioReco.voltageToEfieldConverterPerChannelGroup') + self.antenna_provider = None + self.__counter = 0 + self.begin() + + def begin(self, use_MC_direction=False): + """ + Initializes the module + + Parameters + ---------- + use_MC_direction: bool + If True, the MC direction is used for the reconstruction. If False, the reconstructed angles are used. + """ + self.antenna_provider = antennapattern.AntennaPatternProvider() + self.__use_MC_direction = use_MC_direction + + @register_run() + def run(self, evt, station, det): + """ + Performs computation for voltage trace to electric field per channel + + Will provide a deconvoluted (electric field) trace for each channel from the stations input voltage traces + + Parameters + ---------- + evt: event data structure + the event data structure + station: station data structure + the station data structure + det: detector object + the detector object + """ + if self.__use_MC_direction: + if station.get_sim_station() is not None and station.get_sim_station().has_parameter(stnp.zenith): + zenith = station.get_sim_station()[stnp.zenith] + azimuth = station.get_sim_station()[stnp.azimuth] + else: + self.logger.error(f"MC direction requested but no simulation present in station {station.get_id()}") + raise ValueError("MC direction requested but no simulation present") + else: + self.logger.debug("Using reconstructed angles as no simulation present") + zenith = station[stnp.zenith] + azimuth = station[stnp.azimuth] + + use_channels = det.get_channel_ids(station.get_id()) + frequencies = station.get_channel(use_channels[0]).get_frequencies() # assuming that all channels have the same sampling rate and length + + sampling_rate = station.get_channel(use_channels[0]).get_sampling_rate() + + group_ids = select_channels_per_station(det, station.get_id(), station.get_channel_ids()) + for gid, use_channels in group_ids.items(): + pos = [] + for channel_id in use_channels: + pos.append(det.get_relative_position(station.get_id(), channel_id)) + pos = np.array(pos) + pos = np.average(pos, axis=0) + efield_antenna_factor = trace_utilities.get_efield_antenna_factor(station, frequencies, use_channels, det, + zenith, azimuth, self.antenna_provider) + V = np.zeros((len(use_channels), len(frequencies)), dtype=complex) + for i_ch, channel_id in enumerate(use_channels): + V[i_ch] = station.get_channel(channel_id).get_frequency_spectrum() + denom = (efield_antenna_factor[0][0] * efield_antenna_factor[1][1] - efield_antenna_factor[0][1] * efield_antenna_factor[1][0]) + mask = np.abs(denom) != 0 + # solving for electric field using just two orthorgonal antennas + E1 = np.zeros_like(V[0], dtype=complex) + E2 = np.zeros_like(V[0], dtype=complex) + E1[mask] = (V[0] * efield_antenna_factor[1][1] - V[1] * efield_antenna_factor[0][1])[mask] / denom[mask] + E2[mask] = (V[1] - efield_antenna_factor[1][0] * E1)[mask] / efield_antenna_factor[1][1][mask] + denom = (efield_antenna_factor[0][0] * efield_antenna_factor[-1][1] - efield_antenna_factor[0][1] * efield_antenna_factor[-1][0]) + mask = np.abs(denom) != 0 + E1[mask] = (V[0] * efield_antenna_factor[-1][1] - V[-1] * efield_antenna_factor[0][1])[mask] / denom[mask] + E2[mask] = (V[-1] - efield_antenna_factor[-1][0] * E1)[mask] / efield_antenna_factor[-1][1][mask] + + efield = ef.ElectricField(use_channels) + efield.set_position(pos) + efield.set_frequency_spectrum(np.array([np.zeros_like(E1), E1, E2]), sampling_rate) + + efield.set_trace_start_time(station.get_channel(use_channels[0]).get_trace_start_time()) + efield[efp.zenith] = zenith + efield[efp.azimuth] = azimuth + station.add_electric_field(efield) + + def end(self): + pass diff --git a/NuRadioReco/test/test_examples.sh b/NuRadioReco/test/test_examples.sh index 0fde08ab0..4d7521e46 100755 --- a/NuRadioReco/test/test_examples.sh +++ b/NuRadioReco/test/test_examples.sh @@ -34,7 +34,6 @@ rm output_snr.json #python3 noise_trigger_rate.py --ntries 10 # Full reconstruction cd ../.. -python3 FullReconstruction.py 32 example_data/example_data.hdf5 example_data/arianna_station_32.json python3 read_full_CoREAS_shower.py example_data/example_data.hdf5 python3 SimpleMCReconstruction.py python3 CustomHybridDetector.py diff --git a/NuRadioReco/test/tiny_reconstruction/TinyReconstruction.py b/NuRadioReco/test/tiny_reconstruction/TinyReconstruction.py deleted file mode 100755 index a028e05e5..000000000 --- a/NuRadioReco/test/tiny_reconstruction/TinyReconstruction.py +++ /dev/null @@ -1,176 +0,0 @@ -#!/usr/bin/env python3 -import os -import sys -import datetime -import matplotlib -import matplotlib.pyplot as plt -from NuRadioReco.utilities import units -from NuRadioReco.detector import detector -import NuRadioReco.modules.io.coreas.readCoREAS -import NuRadioReco.modules.io.coreas.simulationSelector -import NuRadioReco.modules.efieldToVoltageConverter -import NuRadioReco.modules.ARIANNA.hardwareResponseIncorporator -import NuRadioReco.modules.channelGenericNoiseAdder -import NuRadioReco.modules.trigger.simpleThreshold -import NuRadioReco.modules.channelBandPassFilter -import NuRadioReco.modules.eventTypeIdentifier -import NuRadioReco.modules.channelStopFilter -import NuRadioReco.modules.channelSignalReconstructor -import NuRadioReco.modules.correlationDirectionFitter -import NuRadioReco.modules.voltageToEfieldConverter -import NuRadioReco.modules.electricFieldSignalReconstructor -import NuRadioReco.modules.electricFieldBandPassFilter -import NuRadioReco.modules.voltageToAnalyticEfieldConverter -import NuRadioReco.modules.channelResampler -import NuRadioReco.modules.electricFieldResampler -import NuRadioReco.modules.io.eventWriter - -# Setup logger for this script -import logging -logger = logging.getLogger('NuRadioReco.TinyReconstruction') - -matplotlib.use('agg') -plt.switch_backend('agg') - - -""" -Here, we show an example reconstruction of CoREAS data. A variety of modules -are being used. Please refer to details in the modules themselves. - -Input parameters (all with a default provided) ---------------------- - -Command line input: - python FullReconstruction.py station_id input_file detector_file templates - -station_id: int - station id to be used, default 32 -input_file: str - CoREAS simulation file, default example data -detector_file: str - path to json detector database, default given -template_path: str - path to signal templates, default given - -""" - -dir_path = os.path.dirname(os.path.realpath(__file__)) # get the directory of this file - -try: - station_id = int(sys.argv[1]) # specify station id - input_file = sys.argv[2] # file with coreas simulations -except: - logger.warning("Usage: python FullReconstruction.py station_id input_file detector templates") - station_id = 32 - input_file = os.path.join(dir_path, "../../examples/example_data/example_event.h5") - logger.warning("Using default station {}".format(32)) - -if(station_id == 32): - triggered_channels = [0, 1, 2, 3] - used_channels_efield = [0, 1, 2, 3] - used_channels_fit = [0, 1, 2, 3] - channel_pairs = ((0, 2), (1, 3)) -else: - logger.warning("Default channels not defined for station_id != 32") - -try: - detector_file = sys.argv[3] - logger.info("Using {0} as detector ".format(detector_file)) -except: - - logger.warning("Using default file for detector") - detector_file = os.path.join(dir_path, "../../examples/example_data/arianna_station_32.json") - -det = detector.Detector(json_filename=detector_file) # detector file -det.update(datetime.datetime(2018, 10, 1)) - -# initialize all modules that are needed for processing -# provide input parameters that are to remain constant during processung -readCoREAS = NuRadioReco.modules.io.coreas.readCoREAS.readCoREAS() -readCoREAS.begin([input_file], station_id, n_cores=10, max_distance=None, seed=0) -simulationSelector = NuRadioReco.modules.io.coreas.simulationSelector.simulationSelector() -simulationSelector.begin() -efieldToVoltageConverter = NuRadioReco.modules.efieldToVoltageConverter.efieldToVoltageConverter() -efieldToVoltageConverter.begin(debug=False) -hardwareResponseIncorporator = NuRadioReco.modules.ARIANNA.hardwareResponseIncorporator.hardwareResponseIncorporator() -channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() -channelGenericNoiseAdder.begin(seed=1) -triggerSimulator = NuRadioReco.modules.trigger.simpleThreshold.triggerSimulator() -triggerSimulator.begin() -channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() -channelBandPassFilter.begin() -eventTypeIdentifier = NuRadioReco.modules.eventTypeIdentifier.eventTypeIdentifier() -channelStopFilter = NuRadioReco.modules.channelStopFilter.channelStopFilter() -channelSignalReconstructor = NuRadioReco.modules.channelSignalReconstructor.channelSignalReconstructor() -channelSignalReconstructor.begin(signal_window_start=20 * units.ns, signal_window_length=80 * units.ns, noise_window_start=150 * units.ns, noise_window_length=200 * units.ns) -correlationDirectionFitter = NuRadioReco.modules.correlationDirectionFitter.correlationDirectionFitter() -voltageToEfieldConverter = NuRadioReco.modules.voltageToEfieldConverter.voltageToEfieldConverter() - -electricFieldSignalReconstructor = NuRadioReco.modules.electricFieldSignalReconstructor.electricFieldSignalReconstructor() -electricFieldSignalReconstructor.begin() - -voltageToAnalyticEfieldConverter = NuRadioReco.modules.voltageToAnalyticEfieldConverter.voltageToAnalyticEfieldConverter() -voltageToAnalyticEfieldConverter.begin() - -electricFieldResampler = NuRadioReco.modules.electricFieldResampler.electricFieldResampler() -electricFieldResampler.begin() -electricFieldBandPassFilter = NuRadioReco.modules.electricFieldBandPassFilter.electricFieldBandPassFilter() - -channelResampler = NuRadioReco.modules.channelResampler.channelResampler() -channelResampler.begin() -eventWriter = NuRadioReco.modules.io.eventWriter.eventWriter() -output_filename = "MC_example_station_{}.nur".format(station_id) -eventWriter.begin(output_filename) - - -event_counter = 0 -# Loop over all events in file as initialized in readCoRREAS and perform analysis -for iE, evt in enumerate(readCoREAS.run(detector=det)): - logger.warning("Processing event number {}".format(event_counter)) - logger.info("processing event {:d} with id {:d}".format(iE, evt.get_id())) - station = evt.get_station(station_id) - - if simulationSelector.run(evt, station.get_sim_station(), det): - - efieldToVoltageConverter.run(evt, station, det) - - hardwareResponseIncorporator.run(evt, station, det, sim_to_data=True) - - channelGenericNoiseAdder.run(evt, station, det, type="rayleigh", amplitude=20 * units.mV) - - triggerSimulator.run(evt, station, det, number_concidences=2, threshold=100 * units.mV) - - if station.get_trigger('default_simple_threshold').has_triggered(): - - channelBandPassFilter.run(evt, station, det, passband=[80 * units.MHz, 500 * units.MHz], filter_type='butter', order=10) - - eventTypeIdentifier.run(evt, station, "forced", 'cosmic_ray') - - channelStopFilter.run(evt, station, det) - - channelBandPassFilter.run(evt, station, det, passband=[60 * units.MHz, 600 * units.MHz], filter_type='rectangular') - channelSignalReconstructor.run(evt, station, det) - - hardwareResponseIncorporator.run(evt, station, det) - - correlationDirectionFitter.run(evt, station, det, n_index=1., channel_pairs=channel_pairs) - - voltageToEfieldConverter.run(evt, station, det, use_channels=used_channels_efield) - - electricFieldBandPassFilter.run(evt, station, det, passband=[80 * units.MHz, 300 * units.MHz]) - - electricFieldSignalReconstructor.run(evt, station, det) - - voltageToAnalyticEfieldConverter.run(evt, station, det, use_channels=used_channels_efield, bandpass=[80 * units.MHz, 500 * units.MHz], use_MC_direction=False) - - channelResampler.run(evt, station, det, sampling_rate=1 * units.GHz) - - electricFieldResampler.run(evt, station, det, sampling_rate=1 * units.GHz) - - eventWriter.run(evt) - - event_counter += 1 - if event_counter > 2: - break -nevents = eventWriter.end() -logger.warning("Finished processing, {} events".format(event_counter)) diff --git a/NuRadioReco/test/tiny_reconstruction/compareToReference.py b/NuRadioReco/test/tiny_reconstruction/compareToReference.py deleted file mode 100644 index 09e80d7f2..000000000 --- a/NuRadioReco/test/tiny_reconstruction/compareToReference.py +++ /dev/null @@ -1,121 +0,0 @@ -#!/usr/bin/env python3 -import json -import NuRadioReco.modules.io.eventReader -import argparse -import numpy as np -from NuRadioReco.framework.parameters import stationParameters as stnp -from NuRadioReco.framework.parameters import channelParameters as chp -from NuRadioReco.framework.parameters import electricFieldParameters as efp - - -parser = argparse.ArgumentParser() -parser.add_argument('filename', help='name of the file that should be compared to reference') -parser.add_argument('reference_file', help='name of the reference file') -parser.add_argument('--station_id', default=32) -parser.add_argument('--create_reference', help='create new reference instead of comparing to current reference', action='store_true') -args = parser.parse_args() - -event_reader = NuRadioReco.modules.io.eventReader.eventReader() -event_reader.begin(args.filename) - -parameter_values = {} - -for event in event_reader.run(): - parameter_values[f"{event.get_id():d}"] = { - 'station_parameters': {}, - 'sim_station_parameters': {}, - 'channel_parameters': {}, - 'electric_field_parameters': {} - } - station = event.get_station(args.station_id) - sim_station = station.get_sim_station() - for param_name in stnp: - if station.has_parameter(param_name): - parameter_values[f"{event.get_id():d}"]['station_parameters'][param_name.name] = station.get_parameter(param_name) - else: - parameter_values[f"{event.get_id():d}"]['station_parameters'][param_name.name] = None - if sim_station.has_parameter(param_name): - parameter_values[f"{event.get_id():d}"]['sim_station_parameters'][param_name.name] = sim_station.get_parameter(param_name) - else: - parameter_values[f"{event.get_id():d}"]['sim_station_parameters'][param_name.name] = None - for channel in station.iter_channels(): - for param_name in chp: - if param_name.name not in parameter_values[f"{event.get_id():d}"]['channel_parameters'].keys(): - parameter_values[f"{event.get_id():d}"]['channel_parameters'][param_name.name] = [] - if channel.has_parameter(param_name): - parameter_values[f"{event.get_id():d}"]['channel_parameters'][param_name.name].append(channel.get_parameter(param_name)) - else: - parameter_values[f"{event.get_id():d}"]['channel_parameters'][param_name.name].append(None) - for electric_field in station.get_electric_fields(): - for param_name in efp: - if param_name.name not in parameter_values[f"{event.get_id():d}"]['electric_field_parameters'].keys(): - parameter_values[f"{event.get_id():d}"]['electric_field_parameters'][param_name.name] = [] - if electric_field.has_parameter(param_name): - value = electric_field.get_parameter(param_name) - if isinstance(value, np.ndarray): - value = list(value) - parameter_values[f"{event.get_id():d}"]['electric_field_parameters'][param_name.name].append(value) - else: - parameter_values[f"{event.get_id():d}"]['electric_field_parameters'][param_name.name].append(None) - - -def assertDeepAlmostEqual(expected, actual, *args, **kwargs): - """ - Assert that two complex structures have almost equal contents. - - Compares lists, dicts and tuples recursively. Checks numeric values - using test_case's :py:meth:`unittest.TestCase.assertAlmostEqual` and - checks all other values with :py:meth:`unittest.TestCase.assertEqual`. - Accepts additional positional and keyword arguments and pass those - intact to assertAlmostEqual() (that's how you specify comparison - precision). - - :param test_case: TestCase object on which we can call all of the basic - 'assert' methods. - :type test_case: :py:class:`unittest.TestCase` object - """ - is_root = '__trace' not in kwargs - trace = kwargs.pop('__trace', 'ROOT') - try: - if isinstance(expected, (int, float, complex)): - np.testing.assert_allclose(actual, expected, *args, **kwargs) - elif isinstance(expected, (list, tuple, np.ndarray)): - np.testing.assert_equal(len(expected), len(actual)) - for index in range(len(expected)): - v1, v2 = expected[index], actual[index] - assertDeepAlmostEqual(v1, v2, - __trace=repr(index), *args, **kwargs) - elif isinstance(expected, dict): - for key in actual: - if key not in expected: - print( - 'Key {} found in reconstruction but not in reference. ' - 'This is fine if you just added it, but you should update ' - 'the reference file to include it in future tests!'.format(key) - ) - else: - assertDeepAlmostEqual( - expected[key], - actual[key], - __trace=repr(key), - *args, - **kwargs - ) - else: - np.testing.assert_equal(actual, expected) - except AssertionError as exc: - exc.__dict__.setdefault('traces', []).append(trace) - if is_root: - trace = ' -> '.join(reversed(exc.traces)) - exc = AssertionError("%s\nTRACE: %s" % (exc, trace)) - raise exc - - -if args.create_reference: - with open(args.reference_file, 'w') as f: - json.dump(parameter_values, f) -else: - with open(args.reference_file, 'r') as f: - parameter_reference = json.load(f) - found_error = False - assertDeepAlmostEqual(parameter_reference, parameter_values, rtol=1e-6) diff --git a/NuRadioReco/test/tiny_reconstruction/reference.json b/NuRadioReco/test/tiny_reconstruction/reference.json deleted file mode 100644 index 968d3584b..000000000 --- a/NuRadioReco/test/tiny_reconstruction/reference.json +++ /dev/null @@ -1 +0,0 @@ -{"0": {"station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": null, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": 0.09219726935932945, "zenith": 0.76585400390625, "azimuth": 3.94125048828125, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": null, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": null, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "sim_station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": 1.58489319246e+18, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": null, "zenith": 0.7853981633974483, "azimuth": 3.957853280977215, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": 1.3918747974426e+18, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": 646.2024663, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "channel_parameters": {"zenith": [null, null, null, null], "azimuth": [null, null, null, null], "maximum_amplitude": [0.0741114148547873, 0.07822416425569524, 0.08019365868532256, 0.09219726935932945], "SNR": [{"integrated_power": 0.0, "peak_2_peak_amplitude": 0.004780092071917627, "peak_amplitude": 0.00485887816779479, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.0029781482057526206, "peak_amplitude": 0.0030899574970230446, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.002963497087732319, "peak_amplitude": 0.0029912274254951183, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.003315202825440897, "peak_amplitude": 0.0035027540130119802, "Seckel_2_noise": 5}], "maximum_amplitude_envelope": [0.0789854151249416, 0.0850534211758587, 0.08035990802734218, 0.09425391321371974], "P2P_amplitude": [0.13819712244096005, 0.15451260952068696, 0.15292399225172826, 0.17650992826729145], "cr_xcorrelations": [null, null, null, null], "nu_xcorrelations": [null, null, null, null], "signal_time": [-165.45597857145674, -169.65597857145679, -182.85597857145672, -181.05597857145676], "noise_rms": [0.009497729224153949, 0.008938569434903116, 0.008083395580028708, 0.009418626437748765], "signal_regions": [null, null, null, null], "noise_regions": [null, null, null, null], "signal_time_offset": [null, null, null, null], "signal_receiving_zenith": [null, null, null, null], "signal_ray_type": [null, null, null, null], "signal_receiving_azimuth": [null, null, null, null]}, "electric_field_parameters": {"ray_path_type": [null, null], "polarization_angle": [1.4654396095601694, 1.4895721505795514], "polarization_angle_expectation": [1.3193420263512765, -1.8222506272385168], "signal_energy_fluence": [[0.0, 0.1499228726724943, 13.406681586865194], [0.0, 0.1138089722425281, 17.17484362157034]], "cr_spectrum_slope": [null, -5.857288373436567], "zenith": [0.76585400390625, 0.76585400390625], "azimuth": [3.94125048828125, 3.94125048828125], "signal_time": [-256.0418188450112, null], "nu_vertex_distance": [null, null], "nu_viewing_angle": [null, null], "max_amp_antenna": [null, null], "max_amp_antenna_envelope": [null, null], "reflection_coefficient_theta": [null, null], "reflection_coefficient_phi": [null, null], "cr_spectrum_quadratic_term": [null, 5.585796124606195], "energy_fluence_ratios": [null, null]}}, "1": {"station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": null, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": 0.34265971445728505, "zenith": 0.79, "azimuth": 3.96, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": null, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": null, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "sim_station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": 1.58489319246e+18, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": null, "zenith": 0.7853981633974483, "azimuth": 3.957853280977215, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": 1.3918747974426e+18, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": 646.2024663, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "channel_parameters": {"zenith": [null, null, null, null], "azimuth": [null, null, null, null], "maximum_amplitude": [0.29720895105426465, 0.34265971445728505, 0.28631950709258547, 0.3355886336552229], "SNR": [{"integrated_power": 0.0, "peak_2_peak_amplitude": 0.006919176246584566, "peak_amplitude": 0.007147260712361315, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.0047795421977155994, "peak_amplitude": 0.00483475221965112, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.007958346991753162, "peak_amplitude": 0.00808818541204719, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.006171263679930979, "peak_amplitude": 0.00634102083568216, "Seckel_2_noise": 5}], "maximum_amplitude_envelope": [0.30136927893948445, 0.3430368419418388, 0.2883085667945778, 0.3430413736731565], "P2P_amplitude": [0.5781561921741288, 0.6591005848543372, 0.5651166372487595, 0.6691483403791307], "cr_xcorrelations": [null, null, null, null], "nu_xcorrelations": [null, null, null, null], "signal_time": [583.9440214285432, 585.5440214285434, 570.5440214285434, 575.1440214285433], "noise_rms": [0.009359269057858015, 0.009003414652873262, 0.008341757854147532, 0.009479378300635246], "signal_regions": [null, null, null, null], "noise_regions": [null, null, null, null], "signal_time_offset": [null, null, null, null], "signal_receiving_zenith": [null, null, null, null], "signal_ray_type": [null, null, null, null], "signal_receiving_azimuth": [null, null, null, null]}, "electric_field_parameters": {"ray_path_type": [null, null], "polarization_angle": [1.4164838522395478, 1.4309598632203113], "polarization_angle_expectation": [1.3259385237455839, -1.8156541298442093], "signal_energy_fluence": [[0.0, 7.388558320124647, 305.36881995591784], [0.0, 6.420446226373949, 324.0685176647595]], "cr_spectrum_slope": [null, -6.704957977720079], "zenith": [0.79, 0.79], "azimuth": [3.96, 3.96], "signal_time": [500.04549368634525, null], "nu_vertex_distance": [null, null], "nu_viewing_angle": [null, null], "max_amp_antenna": [null, null], "max_amp_antenna_envelope": [null, null], "reflection_coefficient_theta": [null, null], "reflection_coefficient_phi": [null, null], "cr_spectrum_quadratic_term": [null, -2.2958549567008015], "energy_fluence_ratios": [null, null]}}, "2": {"station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": null, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": 0.29293556379400426, "zenith": 0.78, "azimuth": 3.95, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": null, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": null, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "sim_station_parameters": {"nu_zenith": null, "nu_azimuth": null, "nu_energy": null, "nu_flavor": null, "ccnc": null, "nu_vertex": null, "inelasticity": null, "triggered": null, "cr_energy": 1.58489319246e+18, "cr_zenith": null, "cr_azimuth": null, "channels_max_amplitude": null, "zenith": 0.7853981633974483, "azimuth": 3.957853280977215, "zenith_cr_templatefit": null, "zenith_nu_templatefit": null, "cr_xcorrelations": null, "nu_xcorrelations": null, "station_time": null, "cr_energy_em": 1.3918747974426e+18, "nu_inttype": null, "chi2_efield_time_direction_fit": null, "ndf_efield_time_direction_fit": null, "cr_xmax": 646.2024663, "vertex_2D_fit": null, "distance_correlations": null, "shower_energy": null, "viewing_angles": null}, "channel_parameters": {"zenith": [null, null, null, null], "azimuth": [null, null, null, null], "maximum_amplitude": [0.2565625777544711, 0.2778654985637024, 0.2672414043783895, 0.29293556379400426], "SNR": [{"integrated_power": 0.0, "peak_2_peak_amplitude": 0.007948759818769018, "peak_amplitude": 0.008209754518559166, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.002601223824751355, "peak_amplitude": 0.0026817225129317146, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.008923287811149969, "peak_amplitude": 0.009007262903406938, "Seckel_2_noise": 5}, {"integrated_power": 0.0, "peak_2_peak_amplitude": 0.0068025654072232655, "peak_amplitude": 0.007054720375806741, "Seckel_2_noise": 5}], "maximum_amplitude_envelope": [0.2618281394747253, 0.28224029243565324, 0.26779467480180347, 0.2932426286212969], "P2P_amplitude": [0.5044222411860875, 0.547608219972491, 0.5318703949237792, 0.5851999240317709], "cr_xcorrelations": [null, null, null, null], "nu_xcorrelations": [null, null, null, null], "signal_time": [-102.65597857145667, -133.45597857145674, -148.2559785714567, -147.45597857145674], "noise_rms": [0.008359098917345579, 0.00848077598024109, 0.008989153845587915, 0.008201090181413125], "signal_regions": [null, null, null, null], "noise_regions": [null, null, null, null], "signal_time_offset": [null, null, null, null], "signal_receiving_zenith": [null, null, null, null], "signal_ray_type": [null, null, null, null], "signal_receiving_azimuth": [null, null, null, null]}, "electric_field_parameters": {"ray_path_type": [null, null], "polarization_angle": [1.4786245054415663, 1.5037562254029073], "polarization_angle_expectation": [1.3232105351636871, -1.818382118426106], "signal_energy_fluence": [[0.0, 2.8633662095356187, 335.1320229092165], [0.0, 1.9443851451542769, 431.3306782776532]], "cr_spectrum_slope": [null, -1.715922037957672], "zenith": [0.78, 0.78], "azimuth": [3.95, 3.95], "signal_time": [-193.88791795047757, null], "nu_vertex_distance": [null, null], "nu_viewing_angle": [null, null], "max_amp_antenna": [null, null], "max_amp_antenna_envelope": [null, null], "reflection_coefficient_theta": [null, null], "reflection_coefficient_phi": [null, null], "cr_spectrum_quadratic_term": [null, 0.5209156977104409], "energy_fluence_ratios": [null, null]}}} \ No newline at end of file diff --git a/NuRadioReco/test/tiny_reconstruction/testTinyReconstruction.sh b/NuRadioReco/test/tiny_reconstruction/testTinyReconstruction.sh deleted file mode 100755 index 0ce1fc0c9..000000000 --- a/NuRadioReco/test/tiny_reconstruction/testTinyReconstruction.sh +++ /dev/null @@ -1,7 +0,0 @@ -set -e -cd NuRadioReco/test/tiny_reconstruction -python3 TinyReconstruction.py -python3 compareToReference.py MC_example_station_32.nur reference.json - -# clean up -rm -v MC_example_station_32.nur \ No newline at end of file diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index 23ea274cd..270494b66 100644 --- a/NuRadioReco/utilities/trace_utilities.py +++ b/NuRadioReco/utilities/trace_utilities.py @@ -106,7 +106,7 @@ def get_electric_field_energy_fluence(electric_field_trace, times, signal_window else: f_signal = np.sum(electric_field_trace[:, signal_window_mask] ** 2, axis=1) dt = times[1] - times[0] - if noise_window_mask is not None: + if noise_window_mask is not None and np.sum(noise_window_mask) > 0: f_noise = np.sum(electric_field_trace[:, noise_window_mask] ** 2, axis=1) f_signal -= f_noise * np.sum(signal_window_mask) / np.sum(noise_window_mask) diff --git a/pyproject.toml b/pyproject.toml index 107d2696c..d23348135 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -56,6 +56,7 @@ crflux = "*" pandas = "*" mattak = {git = "https://github.com/RNO-G/mattak"} runtable = {git = "ssh://git@github.com/RNO-G/rnog-runtable.git"} +cr-pulse-interpolator = {git = "https://github.com/nu-radio/cr-pulse-interpolator"} [tool.poetry.extras] documentation = ["Sphinx", "sphinx-rtd-theme", "numpydoc"] @@ -64,3 +65,4 @@ galacticnoise = ['pygdsm'] ift_reco = ['nifty5', 'pypocketfft'] muon_flux_calc = ['MCEq', 'crflux'] RNO_G_DATA = ["mattak", "runtable", "pandas"] +cr_interpolator = ['cr-pulse-interpolator']