Commit 0dbbddaf authored by Moritz Huetten's avatar Moritz Huetten
Browse files

load pedestal information into monitoring data

parent 89ba99dc
......@@ -12,7 +12,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, EventAndMonDataContainer, TelescopePointingContainer, MonitoringCameraContainer, PedestalContainer
from ctapipe.instrument import TelescopeDescription, SubarrayDescription, OpticsDescription, CameraGeometry
__all__ = ['MAGICEventSource']
......@@ -229,7 +229,10 @@ class MAGICEventSource(EventSource):
counter = 0
# Data container - is initialized once, and data is replaced within it after each yield
data = DataContainer()
if not self.is_mc:
data = EventAndMonDataContainer()
else:
data = DataContainer()
# Telescopes with data:
tels_in_file = ["m1", "m2"]
......@@ -246,6 +249,28 @@ class MAGICEventSource(EventSource):
# Setting the new active run (class MarsRun object)
self.current_run = self._set_active_run(run_number)
# Set monitoring data:
if not self.is_mc:
monitoring_data = self.current_run['data'].monitoring_data
for tel_i, tel_id in enumerate(tels_in_file):
monitoring_camera = MonitoringCameraContainer()
pedestal_info = PedestalContainer()
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'])
monitoring_camera.pedestal = pedestal_info
data.mon.tel[tel_i + 1] = monitoring_camera
# Loop over the events
for event_i in range(self.current_run['data'].n_stereo_events):
# Event and run ids
......@@ -567,6 +592,8 @@ class MarsRun:
files satisfying run_file_mask will be used. Defaults to None.
"""
self.n_camera_pixels = 1039
self.run_file_mask = run_file_mask
# Preparing the lists of M1/2 data files
......@@ -598,10 +625,19 @@ class MarsRun:
if len(run_numbers) > 1:
raise ValueError("Run mask corresponds to more than one run: {}".format(run_numbers))
# Reading the event data
# Reading the data
m1_data = self.load_events(self.m1_file_list, self.is_mc, self.n_camera_pixels)
m2_data = self.load_events(self.m2_file_list, self.is_mc, self.n_camera_pixels)
# Getting the event data
self.event_data = dict()
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)
self.event_data['M1'] = m1_data[0]
self.event_data['M2'] = m2_data[0]
# Getting the monitoring data
self.monitoring_data = dict()
self.monitoring_data['M1'] = m1_data[1]
self.monitoring_data['M2'] = m2_data[1]
# Detecting pedestal events
self.pedestal_ids = self._find_pedestal_events()
......@@ -610,7 +646,6 @@ class MarsRun:
# Detecting mono events
self.mono_ids = self._find_mono_events()
self.n_camera_pixels = 1039
@property
def n_events_m1(self):
......@@ -641,9 +676,9 @@ class MarsRun:
return len(self.pedestal_ids['M2'])
@staticmethod
def load_events(file_list, is_mc):
def load_events(file_list, is_mc, n_camera_pixels):
"""
This method loads events from the pre-defiled file and returns them as a dictionary.
This method loads events and monitoring data from the pre-defiled file and returns them as a dictionary.
Parameters
----------
......@@ -678,7 +713,22 @@ class MarsRun:
# run-wise meta information (same for all events)
mars_meta = dict()
# monitoring information (updated from time to time)
monitoring_data = dict()
monitoring_data['badpixelinfo'] = []
monitoring_data['PedestalMJD'] = scipy.array([])
monitoring_data['PedestalFundamental'] = dict()
monitoring_data['PedestalFundamental']['Mean'] = []
monitoring_data['PedestalFundamental']['Rms'] = []
monitoring_data['PedestalFromExtractor'] = dict()
monitoring_data['PedestalFromExtractor']['Mean'] = []
monitoring_data['PedestalFromExtractor']['Rms'] = []
monitoring_data['PedestalFromExtractorRndm'] = dict()
monitoring_data['PedestalFromExtractorRndm']['Mean'] = []
monitoring_data['PedestalFromExtractorRndm']['Rms'] = []
event_data['file_edges'] = [0]
degrees_per_hour = 15.0
......@@ -721,6 +771,18 @@ class MarsRun:
'MReportDrive.fDec'
]
pedestal_array_list = [
'MTimePedestals.fMjd',
'MTimePedestals.fTime.fMilliSec',
'MTimePedestals.fNanoSec',
'MPedPhotFundamental.fArray.fMean',
'MPedPhotFundamental.fArray.fRms',
'MPedPhotFromExtractor.fArray.fMean',
'MPedPhotFromExtractor.fArray.fRms',
'MPedPhotFromExtractorRndm.fArray.fMean',
'MPedPhotFromExtractorRndm.fArray.fRms'
]
# Info only applicable for MC:
mc_list = [
'MMcEvt.fEnergy',
......@@ -755,6 +817,8 @@ class MarsRun:
# Reading meta information:
mars_meta['is_simulation'] = is_mc
# try to read RunHeaders tree (soft fail if not present)
try:
meta_info = input_file['RunHeaders'].arrays(metainfo_array_list)
......@@ -781,7 +845,7 @@ class MarsRun:
raise ValueError(msg)
# Reading the info only contained in real data
if is_simulation == False:
if not is_simulation:
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
......@@ -790,22 +854,50 @@ class MarsRun:
# 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).astype(bool)
except KeyError:
logger.warning("RunHeaders tree not present in file. Cannot read meta information and assume it is a real data run.")
logger.warning("RunHeaders tree not present in file. Cannot read meta information - will assume it is a real data run.")
badpixelinfo = np.zeros(1183)
is_simulation = False
if is_simulation == False:
# try to read Pedestals tree (soft fail if not present)
if not is_simulation:
try:
pedestal_info = input_file['Pedestals'].arrays(pedestal_array_list)
pedestal_mjd = pedestal_info[b'MTimePedestals.fMjd']
pedestal_millisec = pedestal_info[b'MTimePedestals.fTime.fMilliSec']
pedestal_nanosec = pedestal_info[b'MTimePedestals.fNanoSec']
n_pedestals = len(pedestal_mjd)
pedestal_mjd = pedestal_mjd + (pedestal_millisec / 1e3 + pedestal_nanosec / 1e9) / seconds_per_day
monitoring_data['PedestalMJD'] = scipy.concatenate((monitoring_data['PedestalMJD'], pedestal_mjd))
for quantity in ['Mean', 'Rms']:
for i_pedestal in range(n_pedestals):
monitoring_data['PedestalFundamental'][quantity].append(pedestal_info['MPedPhotFundamental.fArray.f{:s}'.format(quantity).encode()][i_pedestal][:n_camera_pixels])
monitoring_data['PedestalFromExtractor'][quantity].append(pedestal_info['MPedPhotFromExtractor.fArray.f{:s}'.format(quantity).encode()][i_pedestal][:n_camera_pixels])
monitoring_data['PedestalFromExtractorRndm'][quantity].append(pedestal_info['MPedPhotFromExtractorRndm.fArray.f{:s}'.format(quantity).encode()][i_pedestal][:n_camera_pixels])
print(pedestal_info['MPedPhotFromExtractorRndm.fArray.f{:s}'.format('Mean').encode()][i_pedestal])
print(pedestal_info['MPedPhotFundamental.fArray.f{:s}'.format('Mean').encode()][i_pedestal])
except KeyError:
logger.warning("Pedestals tree not present in file.")
# Reading event timing information:
if not is_simulation:
event_times = input_file['Events'].arrays(time_array_list)
# Computing the event arrival time
mjd = event_times[b'MTime.fMjd']
millisec = event_times[b'MTime.fTime.fMilliSec']
nanosec = event_times[b'MTime.fNanoSec']
event_mjd = event_times[b'MTime.fMjd']
event_millisec = event_times[b'MTime.fTime.fMilliSec']
event_nanosec = event_times[b'MTime.fNanoSec']
mjd = mjd + (millisec / 1e3 + nanosec / 1e9) / seconds_per_day
event_mjd = event_mjd + (event_millisec / 1e3 + event_nanosec / 1e9) / seconds_per_day
event_data['MJD'] = scipy.concatenate((event_data['MJD'], event_mjd))
# Reading pointing information (in units of degrees):
if 'MPointingPos.' in input_file['Events']:
......@@ -848,17 +940,17 @@ class MarsRun:
drive_dec_pointing_interpolator = scipy.interpolate.interp1d(drive_mjd, drive_dec, fill_value="extrapolate")
# Interpolating the drive pointing to the event time stamps
pointing_zd = drive_zd_pointing_interpolator(mjd)
pointing_az = drive_az_pointing_interpolator(mjd)
pointing_ra = drive_ra_pointing_interpolator(mjd)
pointing_dec = drive_dec_pointing_interpolator(mjd)
pointing_zd = drive_zd_pointing_interpolator(event_mjd)
pointing_az = drive_az_pointing_interpolator(event_mjd)
pointing_ra = drive_ra_pointing_interpolator(event_mjd)
pointing_dec = drive_dec_pointing_interpolator(event_mjd)
else:
# Not enough data to interpolate the pointing direction.
pointing_zd = scipy.repeat(-1, len(mjd))
pointing_az = scipy.repeat(-1, len(mjd))
pointing_ra = scipy.repeat(-1, len(mjd))
pointing_dec = scipy.repeat(-1, len(mjd))
pointing_zd = scipy.repeat(-1, len(event_mjd))
pointing_az = scipy.repeat(-1, len(event_mjd))
pointing_ra = scipy.repeat(-1, len(event_mjd))
pointing_dec = scipy.repeat(-1, len(event_mjd))
event_data['charge'].append(charge)
event_data['arrival_time'].append(arrival_time)
......@@ -870,9 +962,9 @@ class MarsRun:
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))
if is_simulation == False:
event_data['MJD'] = scipy.concatenate((event_data['MJD'], mjd))
else:
# Reading MC only information:
if is_simulation:
mc_info = input_file['Events'].arrays(mc_list)
# N.B.: For MC, there is only one subrun
event_data['true_energy'] = mc_info[b'MMcEvt.fEnergy']
......@@ -885,7 +977,7 @@ class MarsRun:
event_data['file_edges'].append(len(event_data['trigger_pattern']))
return event_data
return event_data, monitoring_data
def _find_pedestal_events(self):
"""
......
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