Commit 9000e864 authored by Moritz Huetten's avatar Moritz Huetten
Browse files

first working version

parent 0e97a041
Pipeline #65247 canceled with stage
......@@ -65,8 +65,12 @@ class MAGICEventSource(EventSource):
super().__init__(input_url=self.file_list[0], **kwargs)
# Retrieving the list of run numbers corresponding to the data files
run_numbers = list(map(self._get_run_number, self.file_list))
self.run_numbers = np.unique(run_numbers)
run_info = list(map(self._get_run_info_from_name, self.file_list))
run_numbers = [i[0] for i in run_info]
run_types = [i[1] for i in run_info]
self.run_numbers, indices = np.unique(run_numbers, return_index=True)
self.run_types = [run_types[i] for i in indices]
# # Setting up the current run with the first run present in the data
# self.current_run = self._set_active_run(run_number=0)
......@@ -127,10 +131,10 @@ class MAGICEventSource(EventSource):
return is_magic_root_file
@staticmethod
def _get_run_number(file_name):
def _get_run_info_from_name(file_name):
"""
This internal method extracts the run number from
the specified file name.
This internal method extracts the run number and
type (data/MC) from the specified file name.
Parameters
----------
......@@ -143,17 +147,27 @@ class MAGICEventSource(EventSource):
A run number of the file.
"""
mask = r".*\d+_M\d+_(\d+)\.\d+_.*"
parsed_info = re.findall(mask, file_name)
mask_data = r".*\d+_M\d+_(\d+)\.\d+_.*"
mask_mc = r".*_M\d_za\d+to\d+_\d_(\d+)_Y_.*"
mask_mc_alt = r".*_M\d_\d_(\d+)_.*"
if len(re.findall(mask_data, file_name)) > 0:
parsed_info = re.findall(mask_data, file_name)
is_mc = False
elif len(re.findall(mask_mc, file_name)) > 0:
parsed_info = re.findall(mask_mc, file_name)
is_mc = True
else:
parsed_info = re.findall(mask_mc_alt, file_name)
is_mc = True
try:
run_number = int(parsed_info[0])
except IndexError:
raise IndexError('Can not identify the run number of the file {:s}'.format(file_name))
raise IndexError('Can not identify the run number and type (data/MC) of the file {:s}'.format(file_name))
return run_number
return run_number, is_mc
def _set_active_run(self, run_number):
def _set_active_run(self, run_number, is_mc):
"""
This internal method sets the run that will be used for data loading.
......@@ -161,10 +175,12 @@ class MAGICEventSource(EventSource):
----------
run_number: int
The run number to use.
is_mc: Bool
Whether the run is MC or data
Returns
-------
MarsDataRun:
MarsRun:
The run to use
"""
......@@ -174,7 +190,7 @@ class MAGICEventSource(EventSource):
run = dict()
run['number'] = run_number
run['read_events'] = 0
run['data'] = MarsDataRun(run_file_mask=this_run_mask, filter_list=self.file_list)
run['data'] = MarsRun(run_file_mask=this_run_mask, filter_list=self.file_list, is_mc=is_mc)
return run
......@@ -214,15 +230,15 @@ class MAGICEventSource(EventSource):
tels_with_data = {1, 2}
# Loop over the available data runs
for run_number in self.run_numbers:
for index, run_number in enumerate(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']
# Setting the new active run (class MarsDataRun object)
self.current_run = self._set_active_run(run_number)
# Setting the new active run (class MarsRun object)
self.current_run = self._set_active_run(run_number, self.run_types[index])
# Loop over the events
for event_i in range(self.current_run['data'].n_stereo_events):
......@@ -273,9 +289,25 @@ class MAGICEventSource(EventSource):
# data.dl1.tel[i_tel + 1].badpixels = np.array(
# file['dl1/tel' + str(i_tel + 1) + '/badpixels'], dtype=np.bool)
# Adding the event arrival time
time_tmp = Time(event_data['mjd'], scale='utc', format='mjd')
data.trig.gps_time = Time(time_tmp, format='unix', scale='utc', precision=9)
if self.run_types[index] == False:
# Adding the event arrival time
time_tmp = Time(event_data['mjd'], scale='utc', format='mjd')
data.trig.gps_time = Time(time_tmp, format='unix', scale='utc', precision=9)
# Filling weather information
weather = WeatherContainer()
weather.air_temperature = event_data['air_temperature'] * u.deg_C
weather.air_pressure = event_data['air_pressure'] * u.hPa
weather.air_humidity = event_data['air_humidity']
data.weather = weather
else:
data.mc.energy = event_data['true_energy'] * u.GeV
data.mc.alt = (90 - event_data['true_zd']) * u.deg
data.mc.az = event_data['true_az'] * u.deg
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
data.mc.core_x = event_data['true_core_x'] * u.cm
data.mc.core_y = event_data['true_core_y'] * u.cm
# Setting the telescopes with data
data.r0.tels_with_data = tels_with_data
......@@ -283,12 +315,6 @@ class MAGICEventSource(EventSource):
data.dl0.tels_with_data = tels_with_data
data.trig.tels_with_trigger = tels_with_data
# Filling weather information
weather = WeatherContainer()
weather.air_temperature = event_data['air_temperature'] * u.deg_C
weather.air_pressure = event_data['air_pressure'] * u.hPa
weather.air_humidity = event_data['air_humidity']
data.weather = weather
yield data
counter += 1
......@@ -326,7 +352,7 @@ class MAGICEventSource(EventSource):
tels_with_data = {tel_i + 1, }
# Loop over the available data runs
for run_number in self.run_numbers:
for index, run_number in enumerate(self.run_numbers):
# Removing the previously read data run from memory
if self.current_run is not None:
......@@ -334,7 +360,7 @@ class MAGICEventSource(EventSource):
del self.current_run['data']
# Setting the new active run
self.current_run = self._set_active_run(run_number)
self.current_run = self._set_active_run(run_number, self.run_types[index])
if telescope == 'M1':
n_events = self.current_run['data'].n_mono_events_m1
......@@ -388,9 +414,24 @@ class MAGICEventSource(EventSource):
# data.dl1.tel[tel_i + 1].badpixels = np.array(
# file['dl1/tel' + str(i_tel + 1) + '/badpixels'], dtype=np.bool)
# Adding the event arrival time
time_tmp = Time(event_data['mjd'], scale='utc', format='mjd')
data.trig.gps_time = Time(time_tmp, format='unix', scale='utc', precision=9)
if self.run_types[index] == False:
# Adding the event arrival time
time_tmp = Time(event_data['mjd'], scale='utc', format='mjd')
data.trig.gps_time = Time(time_tmp, format='unix', scale='utc', precision=9)
# Filling weather information
weather = WeatherContainer()
weather.air_temperature = event_data['air_temperature'] * u.deg_C
weather.air_pressure = event_data['air_pressure'] * u.hPa
weather.air_humidity = event_data['air_humidity']
data.weather = weather
else:
data.mc.energy = event_data['true_energy'] * u.GeV
data.mc.alt = (90 - event_data['true_zd']) * u.deg
data.mc.az = event_data['true_az'] * u.deg
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
data.mc.core_x = event_data['true_core_x'] * u.cm
data.mc.core_y = event_data['true_core_y'] * u.cm
# Setting the telescopes with data
data.r0.tels_with_data = tels_with_data
......@@ -398,13 +439,6 @@ class MAGICEventSource(EventSource):
data.dl0.tels_with_data = tels_with_data
data.trig.tels_with_trigger = tels_with_data
# Filling weather information
weather = WeatherContainer()
weather.air_temperature = event_data['air_temperature'] * u.deg_C
weather.air_pressure = event_data['air_pressure'] * u.hPa
weather.air_humidity = event_data['air_humidity']
data.weather = weather
yield data
counter += 1
......@@ -441,7 +475,7 @@ class MAGICEventSource(EventSource):
tels_with_data = {tel_i + 1, }
# Loop over the available data runs
for run_number in self.run_numbers:
for index, run_number in enumerate(self.run_numbers):
# Removing the previously read data run from memory
if self.current_run is not None:
......@@ -449,7 +483,7 @@ class MAGICEventSource(EventSource):
del self.current_run['data']
# Setting the new active run
self.current_run = self._set_active_run(run_number)
self.current_run = self._set_active_run(run_number, self.run_types[index])
if telescope == 'M1':
n_events = self.current_run['data'].n_pedestal_events_m1
......@@ -513,12 +547,13 @@ class MAGICEventSource(EventSource):
data.dl0.tels_with_data = tels_with_data
data.trig.tels_with_trigger = tels_with_data
# Filling weather information
weather = WeatherContainer()
weather.air_temperature = event_data['air_temperature'] * u.deg_C
weather.air_pressure = event_data['air_pressure'] * u.hPa
weather.air_humidity = event_data['air_humidity']
data.weather = weather
if self.run_types[index] == False:
# Filling weather information
weather = WeatherContainer()
weather.air_temperature = event_data['air_temperature'] * u.deg_C
weather.air_pressure = event_data['air_pressure'] * u.hPa
weather.air_humidity = event_data['air_humidity']
data.weather = weather
yield data
counter += 1
......@@ -526,314 +561,12 @@ class MAGICEventSource(EventSource):
return
class MAGICEventSourceMC(EventSource):
"""
EventSource for MAGIC calibrated MCs.
"""
_count = 0
def __init__(self, config=None, tool=None, **kwargs):
"""
Constructor
Parameters
----------
config: traitlets.loader.Config
Configuration specified by config file or cmdline arguments.
Used to set traitlet values.
Set to None if no configuration to pass.
tool: ctapipe.core.Tool
Tool executable that is calling this component.
Passes the correct logger to the component.
Set to None if no Tool to pass.
kwargs: dict
Additional parameters to be passed.
NOTE: The file mask of the data to read can be passed with
the 'input_url' parameter.
"""
try:
import uproot
except ImportError:
msg = "The `uproot` python module is required to access the MAGIC data"
self.log.error(msg)
raise
if len(glob.glob(kwargs['input_url'])) > 1:
raise ImportError('MC data can no be loaded with wildcards. Please load them run by run')
self.file_name = kwargs['input_url']
super().__init__(**kwargs)
self.mc_file = MarsMCFile(self.file_name)
# MAGIC telescope positions in m wrt. to the center of CTA simulations
self.magic_tel_positions = {
1: [-27.24, -146.66, 50.00] * u.m,
2: [-96.44, -96.77, 51.00] * u.m
}
# MAGIC telescope description
optics = OpticsDescription.from_name('MAGIC')
geom = CameraGeometry.from_name('MAGICCam')
self.magic_tel_description = TelescopeDescription(name='MAGIC', tel_type='MAGIC', optics=optics, camera=geom)
self.magic_tel_descriptions = {1: self.magic_tel_description, 2: self.magic_tel_description}
self._subarray_info = SubarrayDescription('MAGIC', self.magic_tel_positions, self.magic_tel_descriptions)
@staticmethod
def is_compatible(file_name):
"""
This method checks if the specified file name corresponds
to a MAGIC data file. The result will be True only if all
the files are of ROOT format and contain an 'Events' tree.
Parameters
----------
file_name: str
A file name to check
Returns
-------
bool:
True if the given file is MAGIC data file, False otherwise.
"""
is_magic_root_file = True
try:
import uproot
try:
with uproot.open(file_name) as input_data:
if 'Events' not in input_data:
is_magic_root_file = False
except ValueError:
# uproot raises ValueError if the file is not a ROOT file
is_magic_root_file = False
pass
except ImportError:
if re.match(r'.+_m\d_.+root', file_name.lower()) is None:
is_magic_root_file = False
return is_magic_root_file
@property
def subarray(self):
return self._subarray_info
def _generator(self):
"""
The default event generator. Return the stereo event
generator instance.
Returns
-------
"""
return self._mono_event_generator()
def _mono_event_generator(self):
"""
Mono event generator. Yields DataContainer instances, filled
with the read event data.
Returns
-------
ctapipe.io.containers.DataContainer
"""
counter = 0
# Data container - is initialized once, and data is replaced within it after each yield
data = DataContainer()
data.meta['origin'] = "MAGIC MC"
data.meta['input_url'] = self.input_url
data.meta['is_simulation'] = True
tels_with_data = {self.mc_file.telescope, }
# Loop over the events
for event_i in range(self.mc_file.n_mono_events):
# Event and run ids
event_order_number = self.mc_file.mono_ids[event_i]
event_id = self.mc_file.event_data['daq_event_number'][event_order_number]
obs_id = self.mc_file.run_number
# Reading event data
event_data = self.mc_file.get_mono_event_data(event_i)
# Event counter
data.count = counter
# Setting up the R0 container
data.r0.obs_id = obs_id
data.r0.event_id = event_id
data.r0.tel.clear()
# Setting up the R1 container
data.r1.obs_id = obs_id
data.r1.event_id = event_id
data.r1.tel.clear()
# Setting up the DL0 container
data.dl0.obs_id = obs_id
data.dl0.event_id = event_id
data.dl0.tel.clear()
# Setting up the MC container
data.mc.tel.clear()
# Creating the telescope pointing container
pointing = TelescopePointingContainer()
pointing.azimuth = np.deg2rad(event_data['pointing_az']) * u.rad
pointing.altitude = np.deg2rad(90 - event_data['pointing_zd']) * u.rad
# Adding the pointing container to the event data
data.pointing[self.mc_file.telescope] = pointing
# Adding event charge and peak positions per pixel
data.dl1.tel[self.mc_file.telescope].image = event_data['image']
data.dl1.tel[self.mc_file.telescope].pulse_time = event_data['pulse_time']
# 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.trig.tels_with_trigger = tels_with_data
# mc = data.mc.tel[self.mc_file.telescope]
# mc.dc_to_pe = array_event['laser_calibrations'][tel_id]['calib']
# mc.pedestal = array_event['camera_monitorings'][tel_id]['pedestal']
# mc.reference_pulse_shape = pixel_settings['ref_shape'].astype('float64')
# mc.meta['refstep'] = float(pixel_settings['ref_step'])
# mc.time_slice = float(pixel_settings['time_slice'])
# mc.photo_electron_image = (
# array_event
# .get('photoelectrons', {})
# .get(tel_index, {})
# .get('photoelectrons', np.zeros(n_pixel, dtype='float32'))
# )
data.mc.energy = event_data['true_energy'] * u.GeV
data.mc.alt = (90 - event_data['true_zd']) * u.deg
data.mc.az = event_data['true_az'] * u.deg
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
data.mc.core_x = event_data['true_core_x'] * u.cm
data.mc.core_y = event_data['true_core_y'] * u.cm
yield data
counter += 1
return
def _pedestal_event_generator(self):
"""
Pedestal event generator. Yields DataContainer instances, filled
with the read event data.
Returns
-------
ctapipe.io.containers.DataContainer
"""
counter = 0
# Data container - is initialized once, and data is replaced within it after each yield
data = DataContainer()
data.meta['origin'] = "MAGIC MC"
data.meta['input_url'] = self.input_url
data.meta['is_simulation'] = True
tels_with_data = {self.mc_file.telescope, }
# Loop over the events
for event_i in range(self.mc_file.n_pedestal_events):
# Event and run ids
event_order_number = self.mc_file.pedestal_ids[event_i]
event_id = self.mc_file.event_data['daq_event_number'][event_order_number]
obs_id = self.mc_file.run_number
# Reading event data
event_data = self.mc_file.get_mono_event_data(event_i)
# Event counter
data.count = counter
# Setting up the R0 container
data.r0.obs_id = obs_id
data.r0.event_id = event_id
data.r0.tel.clear()
# Setting up the R1 container
data.r1.obs_id = obs_id
data.r1.event_id = event_id
data.r1.tel.clear()
# Setting up the DL0 container
data.dl0.obs_id = obs_id
data.dl0.event_id = event_id
data.dl0.tel.clear()
# Setting up the MC container
data.mc.tel.clear()
# Creating the telescope pointing container
pointing = TelescopePointingContainer()
pointing.azimuth = np.deg2rad(event_data['pointing_az']) * u.rad
pointing.altitude = np.deg2rad(90 - event_data['pointing_zd']) * u.rad
# Adding the pointing container to the event data
data.pointing[self.mc_file.telescope] = pointing
# Adding event charge and peak positions per pixel
data.dl1.tel[self.mc_file.telescope].image = event_data['image']
data.dl1.tel[self.mc_file.telescope].pulse_time = event_data['pulse_time']
# 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.trig.tels_with_trigger = tels_with_data
# mc = data.mc.tel[self.mc_file.telescope]
# mc.dc_to_pe = array_event['laser_calibrations'][tel_id]['calib']
# mc.pedestal = array_event['camera_monitorings'][tel_id]['pedestal']
# mc.reference_pulse_shape = pixel_settings['ref_shape'].astype('float64')
# mc.meta['refstep'] = float(pixel_settings['ref_step'])
# mc.time_slice = float(pixel_settings['time_slice'])
# mc.photo_electron_image = (
# array_event
# .get('photoelectrons', {})
# .get(tel_index, {})
# .get('photoelectrons', np.zeros(n_pixel, dtype='float32'))
# )
data.mc.energy = event_data['true_energy'] * u.GeV
data.mc.alt = (90 - event_data['true_zd']) * u.deg
data.mc.az = event_data['true_az'] * u.deg
data.mc.shower_primary_id = 1 - event_data['true_shower_primary_id']
data.mc.h_first_int = event_data['true_h_first_int'] * u.m
data.mc.core_x = event_data['true_core_x'] * u.cm
data.mc.core_y = event_data['true_core_y'] * u.cm
yield data
counter += 1
return
class MarsDataRun:
class MarsRun:
"""
This class implements reading of the event data from a single MAGIC data run.
"""
def __init__(self, run_file_mask, filter_list=None):
def __init__(self, run_file_mask, filter_list=None, is_mc=False):
"""
Constructor of the class. Defines the run to use and the camera pixel arrangement.
......@@ -846,9 +579,12 @@ class MarsDataRun:
filter_list: list, optional
A list of files, to which the run_file_mask should be applied. If None, all the
files satisfying run_file_mask will be used. Defaults to None.
is_mc: Bool
Specify whether the run is data or Monte Carlo simulation
"""
self.run_file_mask = run_file_mask
self.is_mc = is_mc
# Preparing the lists of M1/2 data files
file_list = glob.glob(run_file_mask)
......@@ -863,7 +599,12 @@ class MarsDataRun:
self.m2_file_list.sort()
# Retrieving the list of run numbers corresponding to the data files
run_numbers = list(map(self._get_run_number, file_list))
# N.B.: This info reading is redundant and only cross check, if MarsRun
# is not called from within MAGICEventSource. 'run_info' also extracts
# info whether run is MC or data, so is_mc info could also be extracted
# from here.
run_info = list(map(MAGICEventSource._get_run_info_from_name, file_list))
run_numbers = [i[0] for i in run_info]
run_numbers = scipy.unique(run_numbers)
# Checking if a single run is going to be read
......@@ -872,8 +613,8 @@ class MarsDataRun:
# Reading the event data
self.event_data = dict()
self.event_data['M1'] = self.load_events(self.m1_file_list)
self.event_data['M2'] = self.load_events(self.m2_file_list)
self.event_data['M1'] = self.load_events(self.m1_file_list, self.is_mc)
self.event_data['M2'] = self.load_events(self.m2_file_list, self.is_mc)
# Detecting pedestal events
self.pedestal_ids = self._find_pedestal_events()
......@@ -912,32 +653,9 @@ class MarsDataRun:
def n_pedestal_events_m2(self):
return len(self.pedestal_ids['M2'])
@staticmethod
def _get_run_number(file_name):
"""
This internal method extracts the run number from
a specified file name.
Parameters
----------
file_name: str
A file name to process.
Returns
-------
int:
A run number of the file.
"""
mask = r".*\d+_M\d+_(\d+)\.\d+_.*"
parsed_info = re.findall(mask, file_name)