Source code for NuRadioReco.modules.io.coreas.readCoREASStation

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
import logging
logger = logging.getLogger('NuRadioReco.coreas.readCoREASStation')


[docs]class readCoREASStation:
[docs] def begin(self, input_files, station_id, debug=False): """ begin method initialize readCoREAS module Parameters ---------- input_files: input files list of coreas hdf5 files station_id: station id id number of the radio station as defined in detector """ self.__input_files = input_files self.__station_id = station_id self.__current_input_file = 0 self.__current_event = 0 self.__debug = debug
[docs] @register_run() def run(self, detector): """ Reads in all observers in the CoREAS files and returns a new simulated event for each observer with respect to a given detector with a single station. Parameters ---------- detector: Detector object Detector description of the detector that shall be simulated containing one station """ 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() 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 self.__current_input_file += 1
[docs] def end(self): pass