Skip to content
GitLab
Commit 71a37510 authored by Moritz Huetten's avatar Moritz Huetten
Browse files

delete obsolete init_copy.py

parent f36ebcaa
Hide whitespace changes
Inline Side-by-side
ctapipe_io_magic/init_copy.py deleted 100644 → 0
View file @ f36ebcaa
# Event source for MAGIC calibrated data files.
# Requires uproot package (https://github.com/scikit-hep/uproot).
import glob
import re
import scipy
import numpy as np
import scipy.interpolate
from astropy import units as u
from astropy.time import Time
from ctapipe.io.eventsource import EventSource
from ctapipe.io.containers import DataContainer, TelescopePointingContainer, WeatherContainer
from ctapipe.instrument import TelescopeDescription, SubarrayDescription, OpticsDescription, CameraGeometry
__all__ = ['MAGICEventSource']
class MAGICEventSource(EventSource):
"""
EventSource for MAGIC calibrated data.
This class operates with the MAGIC data run-wise. This means that the files
corresponding to the same data run are loaded and processed together.
"""
_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
self.file_list = glob.glob(kwargs['input_url'])
self.file_list.sort()
# EventSource can not handle file wild cards as input_url
# To overcome this we substitute the input_url with first file matching
# the specified file mask.
del kwargs['input_url']
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)
# # Setting up the current run with the first run present in the data
# self.current_run = self._set_active_run(run_number=0)
self.current_run = None
# 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', type='MAGIC', optics=optics, camera=geom)
self.magic_tel_descriptions = {1: self.magic_tel_description, 2: self.magic_tel_description}
self.magic_subarray = SubarrayDescription('MAGIC', self.magic_tel_positions, self.magic_tel_descriptions)
@staticmethod
def is_compatible(file_mask):
"""
This method checks if the specified file mask corresponds
to MAGIC data files. The result will be True only if all
the files are of ROOT format and contain an 'Events' tree.
Parameters
----------
file_mask: str
A file mask to check
Returns
-------
bool:
True if the masked files are MAGIC data runs, False otherwise.
"""
is_magic_root_file = True
file_list = glob.glob(file_mask)
for file_path in file_list:
try:
import uproot
try:
with uproot.open(file_path) 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_path.lower()) is None:
is_magic_root_file = False
return is_magic_root_file
@staticmethod
def _get_run_number(file_name):
"""
This internal method extracts the run number from
the 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)
try:
run_number = int(parsed_info[0])
except IndexError:
raise IndexError('Can not identify the run number of the file {:s}'.format(file_name))
return run_number
def _set_active_run(self, run_number):
"""
This internal method sets the run that will be used for data loading.
Parameters
----------
run_number: int
The run number to use.
Returns
-------
"""
input_path = '/'.join(self.input_url.split('/')[:-1])
this_run_mask = input_path + '/*{:d}*root'.format(run_number)
run = dict()
run['number'] = run_number
run['read_events'] = 0
run['data'] = MarsDataRun(run_file_mask=this_run_mask, filter_list=self.file_list)
return run
def _generator(self):
"""
The default event generator. Return the stereo event
generator instance.
Returns
-------
"""
return self._stereo_event_generator()
def _stereo_event_generator(self):
"""
Stereo event generator. Yields DataContainer instances, filled
with the read event data.
Returns
-------
"""
counter = 0
# Data container - is initialized once, and data is replaced within it after each yield
data = DataContainer()
data.meta['origin'] = "MAGIC"
data.meta['input_url'] = self.input_url
data.meta['is_simulation'] = False
# Telescopes with data:
tels_in_file = ["m1", "m2"]
tels_with_data = {1, 2}
# 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']
# Setting the new active run
self.current_run = self._set_active_run(run_number)
# Loop over the events
for event_i in range(self.current_run['data'].n_stereo_events):
# Event and run ids
event_order_number = self.current_run['data'].stereo_ids[event_i][0]
event_id = self.current_run['data'].event_data['M1']['stereo_event_number'][event_order_number]
obs_id = self.current_run['number']
# Reading event data
event_data = self.current_run['data'].get_stereo_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()
# Filling the DL1 container with the event data
for tel_i, tel_id in enumerate(tels_in_file):
# Creating the telescope pointing container
pointing = TelescopePointingContainer()
pointing.azimuth = np.deg2rad(event_data['{:s}_pointing_az'.format(tel_id)]) * u.rad
pointing.altitude = np.deg2rad(90 - event_data['{:s}_pointing_zd'.format(tel_id)]) * u.rad
pointing.ra = np.deg2rad(event_data['{:s}_pointing_ra'.format(tel_id)]) * u.rad
pointing.dec = np.deg2rad(event_data['{:s}_pointing_dec'.format(tel_id)]) * u.rad
# Adding the pointing container to the event data
data.pointing[tel_i + 1] = pointing
# Adding event charge and peak positions per pixel
data.dl1.tel[tel_i + 1].image = event_data['{:s}_image'.format(tel_id)]
data.dl1.tel[tel_i + 1].pulse_time = event_data['{:s}_pulse_time'.format(tel_id)]
data.dl1.tel[tel_i + 1].badpixels = event_data['{:s}_bad_pixels'.format(tel_id)]
# 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)
# 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
# Setting the instrument sub-array
data.inst.subarray = self.magic_subarray
# Filling weather information
weather = WeatherContainer()
weather.air_temperature = event_data['air_temperature']
weather.air_pressure = event_data['air_pressure']
weather.air_humidity = event_data['air_humidity']
data.weather = weather
yield data
counter += 1
return
def _mono_event_generator(self, telescope):
"""
Mono event generator. Yields DataContainer instances, filled
with the read event data.
Parameters
----------
telescope: str
The telescope for which to return events. Can be either "M1" or "M2".
Returns
-------
"""
counter = 0
telescope = telescope.upper()
# Data container - is initialized once, and data is replaced within it after each yield
data = DataContainer()
data.meta['origin'] = "MAGIC"
data.meta['input_url'] = self.input_url
data.meta['is_simulation'] = False
# Telescopes with data:
tels_in_file = ["M1", "M2"]
if telescope not in tels_in_file:
raise ValueError("Specified telescope {:s} is not in the allowed list {}".format(telescope, tels_in_file))
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']
# Setting the new active run
self.current_run = self._set_active_run(run_number)
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)
# 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()
# 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
pointing.ra = np.deg2rad(event_data['pointing_ra']) * u.rad
pointing.dec = np.deg2rad(event_data['pointing_dec']) * u.rad
# Adding the pointing container to the event data
data.pointing[tel_i + 1] = pointing
# Adding event charge and peak positions per pixel
data.dl1.tel[tel_i + 1].image = event_data['image']
data.dl1.tel[tel_i + 1].pulse_time = event_data['pulse_time']
data.dl1.tel[tel_i + 1].badpixels = event_data['bad_pixels']
# 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)
# 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
# Setting the instrument sub-array
data.inst.subarray = self.magic_subarray
# Filling weather information
weather = WeatherContainer()
weather.air_temperature = event_data['air_temperature']
weather.air_pressure = event_data['air_pressure']
weather.air_humidity = event_data['air_humidity']
data.weather = weather
yield data
counter += 1
return
def _pedestal_event_generator(self, telescope):
"""
Pedestal event generator. Yields DataContainer instances, filled
with the read event data.
Parameters
----------
telescope: str
The telescope for which to return events. Can be either "M1" or "M2".
Returns
-------
"""
counter = 0
telescope = telescope.upper()
# Data container - is initialized once, and data is replaced within it after each yield
data = DataContainer()
data.meta['origin'] = "MAGIC"
data.meta['input_url'] = self.input_url
data.meta['is_simulation'] = False
# Telescopes with data:
tels_in_file = ["M1", "M2"]
if telescope not in tels_in_file:
raise ValueError("Specified telescope {:s} is not in the allowed list {}".format(telescope, tels_in_file))
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']
# Setting the new active run
self.current_run = self._set_active_run(run_number)
if telescope == 'M1':
n_events = self.current_run['data'].n_pedestal_events_m1
else:
n_events = self.current_run['data'].n_pedestal_events_m2
# Loop over the events
for event_i in range(n_events):
# Event and run ids
event_order_number = self.current_run['data'].pedestal_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_pedestal_event_data(event_i, telescope=telescope)
# 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()
# 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
pointing.ra = np.deg2rad(event_data['pointing_ra']) * u.rad
pointing.dec = np.deg2rad(event_data['pointing_dec']) * u.rad
# Adding the pointing container to the event data
data.pointing[tel_i + 1] = pointing
# Adding event charge and peak positions per pixel
data.dl1.tel[tel_i + 1].image = event_data['image']
data.dl1.tel[tel_i + 1].pulse_time = event_data['pulse_time']
data.dl1.tel[tel_i + 1].badpixels = event_data['bad_pixels']
# 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)
# 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
# Setting the instrument sub-array
data.inst.subarray = self.magic_subarray
# Filling weather information
weather = WeatherContainer()
weather.air_temperature = event_data['air_temperature']
weather.air_pressure = event_data['air_pressure']
weather.air_humidity = event_data['air_humidity']
data.weather = weather
yield data
counter += 1
return
class MarsDataRun:
"""
This class implements reading of the event data from a single MAGIC data run.
"""
def __init__(self, run_file_mask, filter_list=None):
"""
Constructor of the class. Defines the run to use and the camera pixel arrangement.
Parameters
----------
run_file_mask: str
A path mask for files belonging to the run. Must correspond to a single run
or an exception will be raised. Must correspond to calibrated ("sorcerer"-level)
data.
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.
"""
self.run_file_mask = run_file_mask
# Preparing the lists of M1/2 data files
file_list = glob.glob(run_file_mask)
# Filtering out extra files if necessary
if filter_list is not None:
file_list = list(set(file_list) & set(filter_list))
self.m1_file_list = list(filter(lambda name: '_M1_' in name, file_list))
self.m2_file_list = list(filter(lambda name: '_M2_' in name, file_list))
self.m1_file_list.sort()
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))
run_numbers = scipy.unique(run_numbers)
# Checking if a single run is going to be read
if len(run_numbers) > 1:
raise ValueError("Run mask corresponds to more than one run: {}".format(run_numbers))
# 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)
# Detecting pedestal events
self.pedestal_ids = self._find_pedestal_events()
# Detecting stereo events
self.stereo_ids = self._find_stereo_events()
# Detecting mono events
self.mono_ids = self._find_mono_events()
self.n_camera_pixels = 1039
@property
def n_events_m1(self):
return len(self.event_data['M1']['MJD'])
@property
def n_events_m2(self):
return len(self.event_data['M2']['MJD'])
@property
def n_stereo_events(self):
return len(self.stereo_ids)
@property
def n_mono_events_m1(self):
return len(self.mono_ids['M1'])
@property
def n_mono_events_m2(self):
return len(self.mono_ids['M2'])
@property
def n_pedestal_events_m1(self):
return len(self.pedestal_ids['M1'])
@property