Commit f36ebcaa authored by Moritz Huetten's avatar Moritz Huetten
Browse files

Merge remote-tracking branch 'origin/moritz_moredatainfo'

Conflicts:
	.gitignore
	ctapipe_io_magic/__init__.py

Includes now:
- MARS meta information (yet incomplete)
- Weather information
- Different pointing information from each telescope
- Bad pixel information (possibly, not yet correct)
parents 404d2cb9 a2646ea6
Pipeline #64949 failed with stage
/.project
/build/
/ctapipe_io_magic.egg-info/
/dist/
/.project
/.pydevproject
/ctapipe_io_magic.egg-info/
ctapipe_io_magic.egg-info/PKG-INFO
ctapipe_io_magic.egg-info/dependency_links.txt
ctapipe_io_magic.egg-info/requires.txt
ctapipe_io_magic.egg-info/SOURCES.txt
ctapipe_io_magic.egg-info/top_level.txt
dist/
ctapipe_io_magic/__pycache__/__init__.cpython-36.pyc
\ No newline at end of file
......@@ -6,6 +6,13 @@ This module implements the *ctapipe* class, needed to read the calibrated data o
Provided *ctapipe* is already installed, the installation can be done like so:
```bash
git clone https://gitlab.mpcdf.mpg.de/ievo/ctapipe_io_magic.git
pip install ./ctapipe_io_magic/
```
This installation via pip (provided, pip is installed) has the advantage to be nicely controlled for belonging to a given conda environment (and to be uninstalled). Alternatively, do
```bash
git clone https://gitlab.mpcdf.mpg.de/ievo/ctapipe_io_magic.git
cd ctapipe_io_magic
......
......@@ -11,7 +11,7 @@ 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
from ctapipe.io.containers import DataContainer, TelescopePointingContainer, WeatherContainer
from ctapipe.instrument import TelescopeDescription, SubarrayDescription, OpticsDescription, CameraGeometry
__all__ = ['MAGICEventSource']
......@@ -205,9 +205,6 @@ class MAGICEventSource(EventSource):
# 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"]
......@@ -233,6 +230,8 @@ class MAGICEventSource(EventSource):
# Reading event data
event_data = self.current_run['data'].get_stereo_event_data(event_i)
data.meta = event_data['mars_meta']
# Event counter
data.count = counter
......@@ -256,10 +255,10 @@ class MAGICEventSource(EventSource):
for tel_i, tel_id in enumerate(tels_in_file):
# 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
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
......@@ -267,6 +266,7 @@ class MAGICEventSource(EventSource):
# 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)
......@@ -279,6 +279,13 @@ class MAGICEventSource(EventSource):
data.r1.tels_with_data = tels_with_data
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
......@@ -305,9 +312,6 @@ class MAGICEventSource(EventSource):
# 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"]
......@@ -343,6 +347,8 @@ class MAGICEventSource(EventSource):
# 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
......@@ -375,6 +381,7 @@ class MAGICEventSource(EventSource):
# 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)
......@@ -388,6 +395,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
yield data
counter += 1
......@@ -413,9 +427,6 @@ class MAGICEventSource(EventSource):
# 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"]
......@@ -451,6 +462,8 @@ class MAGICEventSource(EventSource):
# Reading event data
event_data = self.current_run['data'].get_pedestal_event_data(event_i, telescope=telescope)
data.meta = event_data['mars_meta']
# Event counter
data.count = counter
......@@ -483,6 +496,7 @@ class MAGICEventSource(EventSource):
# 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)
......@@ -495,6 +509,13 @@ class MAGICEventSource(EventSource):
data.r1.tels_with_data = tels_with_data
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
......@@ -945,18 +966,42 @@ class MarsDataRun:
event_data['pointing_ra'] = scipy.array([])
event_data['pointing_dec'] = scipy.array([])
event_data['MJD'] = scipy.array([])
event_data['air_pressure'] = scipy.array([])
event_data['air_humidity'] = scipy.array([])
event_data['air_temperature'] = scipy.array([])
event_data['badpixelinfo'] = []
event_data['mars_meta'] = []
# run-wise meta information (same for all events)
mars_meta = dict()
event_data['file_edges'] = [0]
degrees_per_hour = 15.0
seconds_per_day = 86400.0
seconds_per_hour = 3600.
array_list = ['MCerPhotEvt.fPixels.fPhot', 'MArrivalTime.fData',
'MTriggerPattern.fPrescaled',
'MRawEvtHeader.fStereoEvtNumber', 'MRawEvtHeader.fDAQEvtNumber',
'MTime.fMjd', 'MTime.fTime.fMilliSec', 'MTime.fNanoSec'
]
pointing_array_list = ['MPointingPos.fZd', 'MPointingPos.fAz', 'MPointingPos.fRa', 'MPointingPos.fDec']
pointing_array_list = ['MPointingPos.fZd', 'MPointingPos.fAz', 'MPointingPos.fRa',
'MPointingPos.fDec', 'MPointingPos.fDevZd',
'MPointingPos.fDevAz', 'MPointingPos.fDevHa',
'MPointingPos.fDevDec',
]
drive_array_list = ['MReportDrive.fMjd', 'MReportDrive.fCurrentZd', 'MReportDrive.fCurrentAz',
'MReportDrive.fRa', 'MReportDrive.fDec']
'MReportDrive.fRa', 'MReportDrive.fDec'
]
weather_array_list = ['MTimeWeather.fMjd', 'MTimeWeather.fTime.fMilliSec', 'MTimeWeather.fNanoSec',
'MReportWeather.fPressure', 'MReportWeather.fHumidity', 'MReportWeather.fTemperature']
metainfo_array_list = ['MRawRunHeader.fRunNumber', 'MRawRunHeader.fRunType', 'MRawRunHeader.fSubRunIndex',
'MRawRunHeader.fSourceRA', 'MRawRunHeader.fSourceDEC', 'MRawRunHeader.fTelescopeNumber']
for file_name in file_list:
......@@ -970,21 +1015,85 @@ class MarsDataRun:
trigger_pattern = events[b'MTriggerPattern.fPrescaled']
stereo_event_number = events[b'MRawEvtHeader.fStereoEvtNumber']
# Reading meta information:
meta_info = input_file['RunHeaders'].arrays(metainfo_array_list)
mars_meta['origin'] = "MAGIC"
mars_meta['input_url'] = file_name
mars_meta['number'] = int(meta_info[b'MRawRunHeader.fRunNumber'][0])
#mars_meta['number_subrun'] = int(meta_info[b'MRawRunHeader.fSubRunIndex'][0])
mars_meta['source_ra'] = meta_info[b'MRawRunHeader.fSourceRA'][0] / seconds_per_hour * degrees_per_hour * u.deg
mars_meta['source_dec'] = meta_info[b'MRawRunHeader.fSourceDEC'][0] / seconds_per_hour * u.deg
is_simulation = int(meta_info[b'MRawRunHeader.fRunType'][0])
if is_simulation == 0:
is_simulation = False
elif is_simulation == 256:
is_simulation = True
else:
msg = "Run type (Data or MC) of MAGIC data file not recognised."
self.log.error(msg)
raise
mars_meta['is_simulation'] = is_simulation
# Reading the info only contained in real data
if is_simulation == False:
badpixelinfo = input_file['RunHeaders']['MBadPixelsCam.fArray.fInfo'].array(uproot.asjagged(uproot.asdtype(np.int32))).flatten().reshape((4, 1183), order='F')
# now we have 3 axes:
# 1st axis: Unsuitable pixels
# 2nd axis: Uncalibrated pixels (says why pixel is unsuitable)
# 3rd axis: Bad hardware pixels (says why pixel is unsuitable)
# Each axis cointains a 32bit integer encoding more information about the specific problem, see MARS software, MBADPixelsPix.h
# Here, we however discard this additional information and only grep the "unsuitable" axis.
badpixelinfo = badpixelinfo[1].astype(bool)
else:
badpixelinfo = np.zeros(1183)
# Computing the event arrival time
mjd = events[b'MTime.fMjd']
millisec = events[b'MTime.fTime.fMilliSec']
nanosec = events[b'MTime.fNanoSec']
mjd = mjd + (millisec / 1e3 + nanosec / 1e9) / 86400.0
mjd = mjd + (millisec / 1e3 + nanosec / 1e9) / seconds_per_day
# Reading weather information:
try:
weather_info = input_file['Weather'].arrays(weather_array_list)
weather_time_day = weather_info[b'MTimeWeather.fMjd']
weather_time_millisec = weather_info[b'MTimeWeather.fTime.fMilliSec']
weather_time_nanosec = weather_info[b'MTimeWeather.fNanoSec']
weather_mjd = weather_time_day + (weather_time_millisec/1e3 + weather_time_nanosec/1e9) / seconds_per_day
weather_mjd, weather_indices = np.unique(weather_mjd, return_index = True)
air_pressure_array = weather_info[b'MReportWeather.fPressure'][weather_indices] # hPa
air_humidity_array = weather_info[b'MReportWeather.fHumidity'][weather_indices]
air_temperature_array = weather_info[b'MReportWeather.fTemperature'][weather_indices] # degree celsius
air_pressure_interpolator = scipy.interpolate.interp1d(weather_mjd, air_pressure_array, fill_value="extrapolate")
air_humidity_interpolator = scipy.interpolate.interp1d(weather_mjd, air_humidity_array, fill_value="extrapolate")
air_temperature_interpolator = scipy.interpolate.interp1d(weather_mjd, air_temperature_array, fill_value="extrapolate")
air_pressure = air_pressure_interpolator(mjd) #* u.hPa
air_humidity = air_humidity_interpolator(mjd)
air_temperature = air_temperature_interpolator(mjd) #* u.deg_C
except:
print("Could not find weather information. "
"Set to 0 degree Celsius, 50% humidity, 790hPa ambient pressure.")
air_pressure = scipy.full(len(mjd), 790.) #* u.hPa
air_humidity = scipy.full(len(mjd), 0.5)
air_temperature = scipy.zeros(len(mjd)) #* u.deg_C
# Reading pointing information (in units of degrees):
if 'MPointingPos.' in input_file['Events']:
# Retrieving the telescope pointing direction
pointing = input_file['Events'].arrays(pointing_array_list)
pointing_zd = pointing[b'MPointingPos.fZd']
pointing_az = pointing[b'MPointingPos.fAz']
pointing_ra = pointing[b'MPointingPos.fRa']
pointing_dec = pointing[b'MPointingPos.fDec']
pointing_zd = pointing[b'MPointingPos.fZd'] - pointing[b'MPointingPos.fDevZd']
pointing_az = pointing[b'MPointingPos.fAz'] - pointing[b'MPointingPos.fDevAz']
pointing_ra = (pointing[b'MPointingPos.fRa'] + pointing[b'MPointingPos.fDevHa']) * degrees_per_hour # N.B. the positive sign here, as HA = local sidereal time - ra
pointing_dec = pointing[b'MPointingPos.fDec'] - pointing[b'MPointingPos.fDevDec']
else:
# Getting the telescope drive info
drive = input_file['Drive'].arrays(drive_array_list)
......@@ -992,7 +1101,7 @@ class MarsDataRun:
drive_mjd = drive[b'MReportDrive.fMjd']
drive_zd = drive[b'MReportDrive.fCurrentZd']
drive_az = drive[b'MReportDrive.fCurrentAz']
drive_ra = drive[b'MReportDrive.fRa']
drive_ra = drive[b'MReportDrive.fRa'] * degrees_per_hour
drive_dec = drive[b'MReportDrive.fDec']
# Finding only non-repeating drive entries
......@@ -1031,12 +1140,17 @@ class MarsDataRun:
event_data['charge'].append(charge)
event_data['arrival_time'].append(arrival_time)
event_data['badpixelinfo'].append(badpixelinfo)
event_data['mars_meta'].append(mars_meta)
event_data['trigger_pattern'] = scipy.concatenate((event_data['trigger_pattern'], trigger_pattern))
event_data['stereo_event_number'] = scipy.concatenate((event_data['stereo_event_number'], stereo_event_number))
event_data['stereo_event_number'] = scipy.concatenate((event_data['stereo_event_number'], stereo_event_number)).astype(dtype='int')
event_data['pointing_zd'] = scipy.concatenate((event_data['pointing_zd'], pointing_zd))
event_data['pointing_az'] = scipy.concatenate((event_data['pointing_az'], pointing_az))
event_data['pointing_ra'] = scipy.concatenate((event_data['pointing_ra'], pointing_ra))
event_data['pointing_dec'] = scipy.concatenate((event_data['pointing_dec'], pointing_dec))
event_data['air_pressure'] = scipy.concatenate((event_data['air_pressure'], air_pressure))
event_data['air_humidity'] = scipy.concatenate((event_data['air_humidity'], air_humidity))
event_data['air_temperature'] = scipy.concatenate((event_data['air_temperature'], air_temperature))
event_data['MJD'] = scipy.concatenate((event_data['MJD'], mjd))
......@@ -1242,11 +1356,15 @@ class MarsDataRun:
The output has the following structure:
'image' - photon_content in requested telescope
'pulse_time' - arrival_times in requested telescope
'pointing_az' - pointing azimuth
'pointing_zd' - pointing zenith angle
'pointing_ra' - pointing right ascension
'pointing_dec' - pointing declination
'mjd' - event arrival time
'bad_pixels' - boolean array indicating problematic pixels
'pointing_az' - pointing azimuth [degrees]
'pointing_zd' - pointing zenith angle [degrees]
'pointing_ra' - pointing right ascension [degrees]
'pointing_dec' - pointing declination [degrees]
'mjd' - event arrival time [MJD]
'air_humidity' - relative ambient air humidity
'air_pressure' - ambient air pressure [astropy units]
'air_temperature' - ambient air temperature [astropy units]
"""
file_num = self._get_pedestal_file_num(pedestal_event_num, telescope)
......@@ -1256,15 +1374,21 @@ class MarsDataRun:
photon_content = self.event_data[telescope]['charge'][file_num][id_in_file][:self.n_camera_pixels]
arrival_times = self.event_data[telescope]['arrival_time'][file_num][id_in_file][:self.n_camera_pixels]
bad_pixels = self.event_data[telescope]['badpixelinfo'][file_num][:self.n_camera_pixels]
event_data = dict()
event_data['image'] = photon_content
event_data['pulse_time'] = arrival_times
event_data['bad_pixels'] = bad_pixels
event_data['pointing_az'] = self.event_data[telescope]['pointing_az'][event_id]
event_data['pointing_zd'] = self.event_data[telescope]['pointing_zd'][event_id]
event_data['pointing_ra'] = self.event_data[telescope]['pointing_ra'][event_id]
event_data['pointing_dec'] = self.event_data[telescope]['pointing_dec'][event_id]
event_data['mjd'] = self.event_data[telescope]['MJD'][event_id]
event_data['air_pressure'] = self.event_data[telescope]['air_pressure'][event_id]
event_data['air_humidity'] = self.event_data[telescope]['air_humidity'][event_id]
event_data['air_temperature'] = self.event_data[telescope]['air_temperature'][event_id]
event_data['mars_meta'] = self.event_data[telescope]['mars_meta'][file_num]
return event_data
......@@ -1286,13 +1410,22 @@ class MarsDataRun:
The output has the following structure:
'm1_image' - M1 photon_content
'm1_pulse_time' - M1 arrival_times
'm1_bad_pixels' - boolean array indicating problematic M1 pixels
'm2_image' - M2 photon_content
'm2_pulse_time' - M2 arrival_times
'pointing_az' - pointing azimuth
'pointing_zd' - pointing zenith angle
'pointing_ra' - pointing right ascension
'pointing_dec' - pointing declination
'mjd' - event arrival time
'm2_peak_pos' - M2 arrival_times
'm2_bad_pixels' - boolean array indicating problematic M2 pixels
'm1_pointing_az' - M1 pointing azimuth [degrees]
'm1_pointing_zd' - M1 pointing zenith angle [degrees]
'm1_pointing_ra' - M1 pointing right ascension [degrees]
'm1_pointing_dec' - M1 pointing declination [degrees]
'm2_pointing_az' - M2 pointing azimuth [degrees]
'm2_pointing_zd' - M2 pointing zenith angle [degrees]
'm2_pointing_ra' - M2 pointing right ascension [degrees]
'm2_pointing_dec' - M2 pointing declination [degrees]
'mjd' - event arrival time [MJD]
'air_humidity' - relative ambient air humidity
'air_pressure' - ambient air pressure [astropy units]
'air_temperature' - ambient air temperature [astropy units]
"""
m1_file_num, m2_file_num = self._get_stereo_file_num(stereo_event_num)
......@@ -1304,20 +1437,33 @@ class MarsDataRun:
m1_photon_content = self.event_data['M1']['charge'][m1_file_num][m1_id_in_file][:self.n_camera_pixels]
m1_arrival_times = self.event_data['M1']['arrival_time'][m1_file_num][m1_id_in_file][:self.n_camera_pixels]
m1_bad_pixels = self.event_data['M1']['badpixelinfo'][m1_file_num][:self.n_camera_pixels]
m2_photon_content = self.event_data['M2']['charge'][m2_file_num][m2_id_in_file][:self.n_camera_pixels]
m2_arrival_times = self.event_data['M2']['arrival_time'][m2_file_num][m2_id_in_file][:self.n_camera_pixels]
m2_bad_pixels = self.event_data['M2']['badpixelinfo'][m2_file_num][:self.n_camera_pixels]
event_data = dict()
event_data['m1_image'] = m1_photon_content
event_data['m1_pulse_time'] = m1_arrival_times
event_data['m1_bad_pixels'] = m1_bad_pixels
event_data['m2_image'] = m2_photon_content
event_data['m2_pulse_time'] = m2_arrival_times
event_data['pointing_az'] = self.event_data['M1']['pointing_az'][m1_id]
event_data['pointing_zd'] = self.event_data['M1']['pointing_zd'][m1_id]
event_data['pointing_ra'] = self.event_data['M1']['pointing_ra'][m1_id]
event_data['pointing_dec'] = self.event_data['M1']['pointing_dec'][m1_id]
event_data['m2_bad_pixels'] = m2_bad_pixels
event_data['m1_pointing_az'] = self.event_data['M1']['pointing_az'][m1_id]
event_data['m1_pointing_zd'] = self.event_data['M1']['pointing_zd'][m1_id]
event_data['m1_pointing_ra'] = self.event_data['M1']['pointing_ra'][m1_id]
event_data['m1_pointing_dec'] = self.event_data['M1']['pointing_dec'][m1_id]
event_data['m2_pointing_az'] = self.event_data['M2']['pointing_az'][m2_id]
event_data['m2_pointing_zd'] = self.event_data['M2']['pointing_zd'][m2_id]
event_data['m2_pointing_ra'] = self.event_data['M2']['pointing_ra'][m2_id]
event_data['m2_pointing_dec'] = self.event_data['M2']['pointing_dec'][m2_id]
# get information identical for both telescopes from M1:
event_data['mjd'] = self.event_data['M1']['MJD'][m1_id]
event_data['air_pressure'] = self.event_data['M1']['air_pressure'][m1_id]
event_data['air_humidity'] = self.event_data['M1']['air_humidity'][m1_id]
event_data['air_temperature'] = self.event_data['M1']['air_temperature'][m1_id]
event_data['mars_meta'] = self.event_data['M1']['mars_meta'][m1_file_num]
return event_data
......@@ -1342,11 +1488,15 @@ class MarsDataRun:
The output has the following structure:
'image' - photon_content in requested telescope
'pulse_time' - arrival_times in requested telescope
'pointing_az' - pointing azimuth
'pointing_zd' - pointing zenith angle
'pointing_ra' - pointing right ascension
'pointing_dec' - pointing declination
'mjd' - event arrival time
'bad_pixels' - boolean array indicating problematic pixels
'pointing_az' - pointing azimuth [degrees]
'pointing_zd' - pointing zenith angle [degrees]
'pointing_ra' - pointing right ascension [degrees]
'pointing_dec' - pointing declination [degrees]
'mjd' - event arrival time [MJD]
'air_humidity' - relative ambient air humidity
'air_pressure' - ambient air pressure [astropy units]
'air_temperature' - ambient air temperature [astropy units]
"""
file_num = self._get_mono_file_num(mono_event_num, telescope)
......@@ -1356,15 +1506,21 @@ class MarsDataRun:
photon_content = self.event_data[telescope]['charge'][file_num][id_in_file][:self.n_camera_pixels]
arrival_times = self.event_data[telescope]['arrival_time'][file_num][id_in_file][:self.n_camera_pixels]
bad_pixels = self.event_data[telescope]['badpixelinfo'][file_num][:self.n_camera_pixels]
event_data = dict()
event_data['image'] = photon_content
event_data['pulse_time'] = arrival_times
event_data['bad_pixels'] = bad_pixels
event_data['pointing_az'] = self.event_data[telescope]['pointing_az'][event_id]
event_data['pointing_zd'] = self.event_data[telescope]['pointing_zd'][event_id]
event_data['pointing_ra'] = self.event_data[telescope]['pointing_ra'][event_id]
event_data['pointing_dec'] = self.event_data[telescope]['pointing_dec'][event_id]
event_data['mjd'] = self.event_data[telescope]['MJD'][event_id]
event_data['air_pressure'] = self.event_data[telescope]['air_pressure'][event_id]
event_data['air_humidity'] = self.event_data[telescope]['air_humidity'][event_id]
event_data['air_temperature'] = self.event_data[telescope]['air_temperature'][event_id]
event_data['mars_meta'] = self.event_data[telescope]['mars_meta'][file_num]
return event_data
......
# 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