Skip to content

Commit

Permalink
change corsika to corsika_evt to avoid confusion with open hdf5 file,…
Browse files Browse the repository at this point in the history
… improve docstring
  • Loading branch information
lpyras committed Jun 3, 2024
1 parent 2040847 commit 4e807d1
Show file tree
Hide file tree
Showing 2 changed files with 65 additions and 26 deletions.
31 changes: 22 additions & 9 deletions NuRadioReco/modules/io/coreas/coreas.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
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 cr_pulse_interpolator.signal_interpolation_fourier
import cr_pulse_interpolator.self.__corsika_evt
import cr_pulse_interpolator.interpolation_fourier
import logging
import copy
Expand Down Expand Up @@ -308,7 +308,7 @@ def make_sim_shower(corsika, observer=None, detector=None, station_id=None):
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
Expand Down Expand Up @@ -495,11 +495,11 @@ class coreasInterpolator:
Parameters
----------
event : NuRadioReco event object containing the CoREAS output (use read_CORSIKA7() to read the file)
corsika_evt : NuRadio event object
use read_CORSIKA7() to create the event object containing the CoREAS output
"""

def __init__(self, evt):
def __init__(self, corsika_evt):
self.electric_field_on_sky = None
self.efield_times = None
self.obs_positions_geo = None
Expand All @@ -509,9 +509,9 @@ def __init__(self, evt):
self.star_radius = None
self.geo_star_radius = None

self.evt = evt
self.sim_station = evt.get_station(0).get_sim_station()
self.shower = evt.get_first_sim_shower() # there should only be one simulated shower
self.corsika_evt = corsika_evt
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.star_shape_initilized = False
self.efield_interpolator_initilized = False
self.fluence_interpolator_initilized = False
Expand Down Expand Up @@ -627,7 +627,7 @@ def initialize_efield_interpolator(self, interp_lowfreq, interp_highfreq):
self.efield_interpolator = -1
else:
logger.info(f'electric field interpolation with lowfreq {interp_lowfreq/units.MHz} MHz and highfreq {interp_highfreq/units.MHz} MHz')
self.efield_interpolator = cr_pulse_interpolator.signal_interpolation_fourier.interp2d_signal(
self.efield_interpolator = cr_pulse_interpolator.self.__corsika_evt.interp2d_signal(
self.obs_positions_vxB_vxvxB[:, 0],
self.obs_positions_vxB_vxvxB[:, 1],
self.electric_field_on_sky,
Expand Down Expand Up @@ -771,6 +771,19 @@ def get_interp_efield_value(self, position_on_ground, core):
def get_interp_fluence_value(self, position_on_ground, core):
"""
Accesses the interpolated fluence for a given position on ground
Parameters
----------
position_on_ground : np.array (3)
position of the antenna on ground
core : np.array (3)
position of the core on ground
Returns
-------
fluence_interp : float
interpolated fluence value
"""
antenna_position = position_on_ground
z_plane = core[2]
Expand Down
60 changes: 43 additions & 17 deletions NuRadioReco/modules/io/coreas/readCoREASDetector.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
import h5py
import numpy as np
from radiotools import coordinatesystems as cstrafo
import cr_pulse_interpolator.signal_interpolation_fourier
import cr_pulse_interpolator.interpolation_fourier
import matplotlib.pyplot as plt
from NuRadioReco.modules.base.module import register_run
Expand All @@ -25,6 +24,29 @@
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
Expand All @@ -35,7 +57,7 @@ def get_random_core_positions(xmin, xmax, ymin, ymax, n_cores, seed=None):

def apply_hanning(efields):
"""
Apply a half hann window to the electric field in the time domain
Apply a half hann window to the electric field in the time domain. This smoothens the edges to avoid ringing effects.
Parameters
----------
Expand Down Expand Up @@ -68,6 +90,11 @@ def select_channels_per_station(det, station_id, requested_channel_ids):
The station id to select channels from
requested_channel_ids : list
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:
Expand All @@ -90,36 +117,35 @@ def __init__(self):
self.__t = 0
self.__t_event_structure = 0
self.__t_per_event = 0
self.__corsika = None
self.__corsika_evt = None
self.__interp_lowfreq = None
self.__interp_highfreq = 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, debug=False):
def begin(self, input_file, interp_lowfreq=30*units.MHz, interp_highfreq=1000*units.MHz, log_level=logging.INFO):
"""
begin method
initialize readCoREAS module
begin method, initialize readCoREASDetector module
Parameters
----------
input_files: input files
list of coreas hdf5 files
interp_lowfreq: float (default = 30)
input_files: str
coreas hdf5 file
interp_lowfreq: float (default = 30 MHz)
lower frequency for the bandpass filter in interpolation, should be broader than the sensetivity band of the detector
interp_highfreq: float (default = 1000)
interp_highfreq: float (default = 1000 MHz)
higher frequency for the bandpass filter in interpolation, should be broader than the sensetivity band of the detector
"""
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 = coreas.read_CORSIKA7(input_file)
self.logger.info(f"using coreas simulation {input_file} with E={self.__corsika.get_first_sim_shower().get_parameter(shp.energy):.2g}eV, zenith angle = {self.__corsika.get_first_sim_shower().get_parameter(shp.zenith) / units.deg:.0f}deg")
self.__corsika_evt = coreas.read_CORSIKA7(input_file)
self.logger.info(f"using coreas simulation {input_file} with E={self.__corsika_evt.get_first_sim_shower().get_parameter(shp.energy):.2g}eV, 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.debug = debug
self.coreasInterpolator = coreas.coreasInterpolator(self.__corsika)
self.coreasInterpolator = coreas.coreasInterpolator(self.__corsika_evt)
self.coreasInterpolator.initialize_efield_interpolator(self.__interp_lowfreq, self.__interp_highfreq)

@register_run()
Expand All @@ -144,13 +170,13 @@ def run(self, detector, core_position_list=[], selected_station_ids=[], selected

for iCore, core in enumerate(core_position_list):
t = time.time()
evt = NuRadioReco.framework.event.Event(self.__corsika.get_run_number(), iCore) # create empty event
sim_shower = self.__corsika.get_first_sim_shower()
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.get_station(0).get_sim_station()
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)
Expand Down

0 comments on commit 4e807d1

Please sign in to comment.