Commit 99650167 authored by Alessio Berti's avatar Alessio Berti
Browse files

Select single run number.

parent 1d72e17d
Pipeline #116184 failed with stage
in 7 seconds
......@@ -478,162 +478,159 @@ class MAGICEventSource(EventSource):
tel_i = tels_in_file.index(telescope)
tels_with_data = {tel_i + 1, }
# Loop over the available data runs
for run_number in self.run_numbers:
# Removing the previously read data run from memory
if self.current_run is not None:
if 'data' in self.current_run:
del self.current_run['data']
# Removing the previously read data run from memory
if self.current_run is not None:
if 'data' in self.current_run:
del self.current_run['data']
# Setting the new active run
self.current_run = self._set_active_run(run_number)
# Setting the new active run
self.current_run = self._set_active_run(self.run_numbers)
# Set monitoring data:
if not self.is_mc:
monitoring_data = self.current_run['data'].monitoring_data
monitoring_camera = MonitoringCameraContainer()
pedestal_info = PedestalContainer()
badpixel_info = PixelStatusContainer()
time_tmp = Time(monitoring_data['M{:d}'.format(
tel_i + 1)]['PedestalMJD'], scale='utc', format='mjd')
pedestal_info.sample_time = Time(
time_tmp, format='unix', scale='utc', precision=9)
pedestal_info.n_events = 500 # hardcoded number of pedestal events averaged over
pedestal_info.charge_mean = []
pedestal_info.charge_mean.append(
monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Mean'])
pedestal_info.charge_mean.append(
monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Mean'])
pedestal_info.charge_mean.append(monitoring_data['M{:d}'.format(
tel_i + 1)]['PedestalFromExtractorRndm']['Mean'])
pedestal_info.charge_std = []
pedestal_info.charge_std.append(
monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Rms'])
pedestal_info.charge_std.append(
monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Rms'])
pedestal_info.charge_std.append(
monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractorRndm']['Rms'])
t_range = Time(monitoring_data['M{:d}'.format(
tel_i + 1)]['badpixelinfoMJDrange'], scale='utc', format='mjd')
badpixel_info.hardware_failing_pixels = monitoring_data['M{:d}'.format(
tel_i + 1)]['badpixelinfo']
badpixel_info.sample_time_range = t_range
monitoring_camera.pedestal = pedestal_info
monitoring_camera.pixel_status = badpixel_info
data.mon.tels_with_data = tels_with_data
data.mon.tel[tel_i + 1] = monitoring_camera
#fdp (mono version not fully tested)
else:
data.mcheader.num_showers = self.current_run['data'].mcheader_data[telescope]['sim_nevents'] # total, including reuse
data.mcheader.shower_reuse = self.current_run['data'].mcheader_data[telescope]['sim_reuse']
data.mcheader.energy_range_min = (self.current_run['data'].mcheader_data[telescope]['sim_emin']).to(u.TeV) # GeV->TeV
data.mcheader.energy_range_max = (self.current_run['data'].mcheader_data[telescope]['sim_emax']).to(u.TeV) # GeV->TeV
data.mcheader.spectral_index = self.current_run['data'].mcheader_data[telescope]['sim_eslope']
data.mcheader.max_scatter_range = (self.current_run['data'].mcheader_data[telescope]['sim_max_impact']).to(u.m) # cm->m
data.mcheader.max_viewcone_radius = (self.current_run['data'].mcheader_data[telescope]['sim_conesemiangle']).to(u.deg) # deg->deg
if telescope == 'M1':
n_events = self.current_run['data'].n_mono_events_m1
else:
n_events = self.current_run['data'].n_mono_events_m2
# Loop over the events
for event_i in range(n_events):
# Event and run ids
event_order_number = self.current_run['data'].mono_ids[telescope][event_i]
event_id = self.current_run['data'].event_data[telescope]['stereo_event_number'][event_order_number]
obs_id = self.current_run['number']
# Reading event data
event_data = self.current_run['data'].get_mono_event_data(
event_i, telescope=telescope)
data.meta = event_data['mars_meta']
# Set monitoring data:
if not self.is_mc:
# Event counter
data.count = counter
data.index.obs_id = obs_id
data.index.event_id = event_id
monitoring_data = self.current_run['data'].monitoring_data
# Setting up the R0 container
data.r0.tel.clear()
data.r0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
monitoring_camera = MonitoringCameraContainer()
pedestal_info = PedestalContainer()
badpixel_info = PixelStatusContainer()
# Setting up the R1 container
data.r1.tel.clear()
data.r1.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
time_tmp = Time(monitoring_data['M{:d}'.format(
tel_i + 1)]['PedestalMJD'], scale='utc', format='mjd')
pedestal_info.sample_time = Time(
time_tmp, format='unix', scale='utc', precision=9)
pedestal_info.n_events = 500 # hardcoded number of pedestal events averaged over
pedestal_info.charge_mean = []
pedestal_info.charge_mean.append(
monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Mean'])
pedestal_info.charge_mean.append(
monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Mean'])
pedestal_info.charge_mean.append(monitoring_data['M{:d}'.format(
tel_i + 1)]['PedestalFromExtractorRndm']['Mean'])
pedestal_info.charge_std = []
pedestal_info.charge_std.append(
monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFundamental']['Rms'])
pedestal_info.charge_std.append(
monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractor']['Rms'])
pedestal_info.charge_std.append(
monitoring_data['M{:d}'.format(tel_i + 1)]['PedestalFromExtractorRndm']['Rms'])
# Setting up the DL0 container
data.dl0.tel.clear()
data.dl0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
t_range = Time(monitoring_data['M{:d}'.format(
tel_i + 1)]['badpixelinfoMJDrange'], scale='utc', format='mjd')
# Setting up the DL1 container
data.dl1.tel.clear()
badpixel_info.hardware_failing_pixels = monitoring_data['M{:d}'.format(
tel_i + 1)]['badpixelinfo']
badpixel_info.sample_time_range = t_range
# Creating the telescope pointing container
pointing = PointingContainer()
pointing_tel = TelescopePointingContainer()
monitoring_camera.pedestal = pedestal_info
monitoring_camera.pixel_status = badpixel_info
pointing_tel.azimuth = np.deg2rad(
event_data['pointing_az']) * u.rad
pointing_tel.altitude = np.deg2rad(
90 - event_data['pointing_zd']) * u.rad
#pointing.ra = np.deg2rad(event_data['pointing_ra']) * u.rad
#pointing.dec = np.deg2rad(event_data['pointing_dec']) * u.rad
pointing.tel[tel_i + 1] = pointing_tel
pointing.array_azimuth = np.deg2rad(event_data['pointing_az']) * u.rad
pointing.array_altitude = np.deg2rad(90 - event_data['pointing_zd']) * u.rad
pointing.array_ra = np.deg2rad(event_data['pointing_ra']) * u.rad
pointing.array_dec = np.deg2rad(90 - event_data['pointing_dec']) * u.rad
data.pointing = pointing
data.mon.tels_with_data = tels_with_data
data.mon.tel[tel_i + 1] = monitoring_camera
# Adding event charge and peak positions per pixel
data.dl1.tel[tel_i + 1].image = event_data['image']
data.dl1.tel[tel_i + 1].peak_time = event_data['pulse_time']
#fdp (mono version not fully tested)
else:
data.mcheader.num_showers = self.current_run['data'].mcheader_data[telescope]['sim_nevents'] # total, including reuse
data.mcheader.shower_reuse = self.current_run['data'].mcheader_data[telescope]['sim_reuse']
data.mcheader.energy_range_min = (self.current_run['data'].mcheader_data[telescope]['sim_emin']).to(u.TeV) # GeV->TeV
data.mcheader.energy_range_max = (self.current_run['data'].mcheader_data[telescope]['sim_emax']).to(u.TeV) # GeV->TeV
data.mcheader.spectral_index = self.current_run['data'].mcheader_data[telescope]['sim_eslope']
data.mcheader.max_scatter_range = (self.current_run['data'].mcheader_data[telescope]['sim_max_impact']).to(u.m) # cm->m
data.mcheader.max_viewcone_radius = (self.current_run['data'].mcheader_data[telescope]['sim_conesemiangle']).to(u.deg) # deg->deg
if not self.is_mc:
# Adding the event arrival time
time_tmp = Time(
event_data['mjd'], scale='utc', format='mjd')
data.trigger.time = Time(
time_tmp, format='unix', scale='utc', precision=9)
else:
data.mc.energy = event_data['true_energy'] * u.GeV
data.mc.alt = (np.pi/2 - event_data['true_zd']) * u.rad
# check meaning of 7deg transformation (I.Vovk)
data.mc.az = -1 * \
(event_data['true_az'] - np.deg2rad(180 - 7)) * u.rad
data.mc.shower_primary_id = 1 - \
event_data['true_shower_primary_id']
data.mc.h_first_int = event_data['true_h_first_int'] * u.cm
# adding a 7deg rotation between the orientation of corsika (x axis = magnetic north) and MARS (x axis = geographical north) frames
# magnetic north is 7 deg westward w.r.t. geographical north
rot_corsika = 7 *u.deg
data.mc.core_x = (event_data['true_core_x']*np.cos(rot_corsika) - event_data['true_core_y']*np.sin(rot_corsika))* u.cm
data.mc.core_y = (event_data['true_core_x']*np.sin(rot_corsika) + event_data['true_core_y']*np.cos(rot_corsika))* u.cm
if telescope == 'M1':
n_events = self.current_run['data'].n_mono_events_m1
else:
n_events = self.current_run['data'].n_mono_events_m2
# Loop over the events
for event_i in range(n_events):
# Event and run ids
event_order_number = self.current_run['data'].mono_ids[telescope][event_i]
event_id = self.current_run['data'].event_data[telescope]['stereo_event_number'][event_order_number]
obs_id = self.current_run['number']
# Reading event data
event_data = self.current_run['data'].get_mono_event_data(
event_i, telescope=telescope)
data.meta = event_data['mars_meta']
# Event counter
data.count = counter
data.index.obs_id = obs_id
data.index.event_id = event_id
# Setting up the R0 container
data.r0.tel.clear()
data.r0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
# Setting up the R1 container
data.r1.tel.clear()
data.r1.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
# Setting up the DL0 container
data.dl0.tel.clear()
data.dl0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
# Setting up the DL1 container
data.dl1.tel.clear()
# Creating the telescope pointing container
pointing = PointingContainer()
pointing_tel = TelescopePointingContainer()
pointing_tel.azimuth = np.deg2rad(
event_data['pointing_az']) * u.rad
pointing_tel.altitude = np.deg2rad(
90 - event_data['pointing_zd']) * u.rad
#pointing.ra = np.deg2rad(event_data['pointing_ra']) * u.rad
#pointing.dec = np.deg2rad(event_data['pointing_dec']) * u.rad
pointing.tel[tel_i + 1] = pointing_tel
pointing.array_azimuth = np.deg2rad(event_data['pointing_az']) * u.rad
pointing.array_altitude = np.deg2rad(90 - event_data['pointing_zd']) * u.rad
pointing.array_ra = np.deg2rad(event_data['pointing_ra']) * u.rad
pointing.array_dec = np.deg2rad(90 - event_data['pointing_dec']) * u.rad
data.pointing = pointing
# Setting the telescopes with data
data.r0.tels_with_data = tels_with_data
data.r1.tels_with_data = tels_with_data
data.dl0.tels_with_data = tels_with_data
data.trigger.tels_with_trigger = tels_with_data
# Adding event charge and peak positions per pixel
data.dl1.tel[tel_i + 1].image = event_data['image']
data.dl1.tel[tel_i + 1].peak_time = event_data['pulse_time']
yield data
counter += 1
if not self.is_mc:
# Adding the event arrival time
time_tmp = Time(
event_data['mjd'], scale='utc', format='mjd')
data.trigger.time = Time(
time_tmp, format='unix', scale='utc', precision=9)
else:
data.mc.energy = event_data['true_energy'] * u.GeV
data.mc.alt = (np.pi/2 - event_data['true_zd']) * u.rad
# check meaning of 7deg transformation (I.Vovk)
data.mc.az = -1 * \
(event_data['true_az'] - np.deg2rad(180 - 7)) * u.rad
data.mc.shower_primary_id = 1 - \
event_data['true_shower_primary_id']
data.mc.h_first_int = event_data['true_h_first_int'] * u.cm
# adding a 7deg rotation between the orientation of corsika (x axis = magnetic north) and MARS (x axis = geographical north) frames
# magnetic north is 7 deg westward w.r.t. geographical north
rot_corsika = 7 *u.deg
data.mc.core_x = (event_data['true_core_x']*np.cos(rot_corsika) - event_data['true_core_y']*np.sin(rot_corsika))* u.cm
data.mc.core_y = (event_data['true_core_x']*np.sin(rot_corsika) + event_data['true_core_y']*np.cos(rot_corsika))* u.cm
# Setting the telescopes with data
data.r0.tels_with_data = tels_with_data
data.r1.tels_with_data = tels_with_data
data.dl0.tels_with_data = tels_with_data
data.trigger.tels_with_trigger = tels_with_data
yield data
counter += 1
return
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment