Commit 910dfbd7 authored by Alessio Berti's avatar Alessio Berti
Browse files

Merge branch 'dev-mhuetten-hotpixels' into 'master'

Dev mhuetten hotpixels

See merge request !9
parents 5a6e2c7a 162644b5
Pipeline #78088 passed with stage
in 10 minutes and 16 seconds
...@@ -70,8 +70,23 @@ pedestal_event_generator = event_source._pedestal_event_generator(telescope='M1' ...@@ -70,8 +70,23 @@ pedestal_event_generator = event_source._pedestal_event_generator(telescope='M1'
... ...
``` ```
#### Features
##### Monitoring data
Monitoring data are saved in `run['data'].monitoring_data` and can also accessed event-wise via the `event.mon` container. Even if they can be accessed event-wise, they are saved only once per run, i.e., identical for all events in a run. If monitoring data is taken several times during a run, the `run['data'].monitoring_data`/`event.mon` sub-containers contain arrays of the quantities taken at the different times together with an array of the time stamps. So far, we have
- Dead pixel information (MARS `RunHeaders.MBadPixelsCam.fArray.fInfo` tree), once per sub-run in `run['data'].monitoring_data['MX']['badpixelinfo']` (with X=1 or X=2) or `event.mon.tel[tel_id].pixel_status`
- Pedestal information from MARS `Pedestals` tree to calculate hot pixels:
- `run['data'].monitoring_data['MX']['PedestalFundamental']`
- `run['data'].monitoring_data['MX']['PedestalFromExtractor']`
- `run['data'].monitoring_data['MX']['PedestalFromExtractorRndm']`
or everyting also in `event.mon.tel[tel_id].pedestal`
Dead pixel and pedestal information are read by `magic-cta-pipe` `MAGIC_Badpixels.py` class.
#### Changelog #### Changelog
- v0.1: Initial version - v0.1: Initial version
- v0.2.0: Unification of data and MC reading - v0.2.0: Unification of data and MC reading
- v0.2.1: Monitoring data (Dead pixel and pedestal information)
\ No newline at end of file
...@@ -12,7 +12,8 @@ import scipy.interpolate ...@@ -12,7 +12,8 @@ import scipy.interpolate
from astropy import units as u from astropy import units as u
from astropy.time import Time from astropy.time import Time
from ctapipe.io.eventsource import EventSource from ctapipe.io.eventsource import EventSource
from ctapipe.io.containers import DataContainer, TelescopePointingContainer from ctapipe.core import Container, Field
from ctapipe.io.containers import DataContainer, EventAndMonDataContainer, TelescopePointingContainer, MonitoringCameraContainer, PedestalContainer
from ctapipe.instrument import TelescopeDescription, SubarrayDescription, OpticsDescription, CameraGeometry from ctapipe.instrument import TelescopeDescription, SubarrayDescription, OpticsDescription, CameraGeometry
__all__ = ['MAGICEventSource'] __all__ = ['MAGICEventSource']
...@@ -20,6 +21,18 @@ __all__ = ['MAGICEventSource'] ...@@ -20,6 +21,18 @@ __all__ = ['MAGICEventSource']
logger = logging.getLogger(__name__) logger = logging.getLogger(__name__)
# MAGIC telescope positions in m wrt. to the center of CTA simulations
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')
magic_tel_description = TelescopeDescription(name='MAGIC', tel_type='MAGIC', optics=optics, camera=geom)
magic_tel_descriptions = {1: magic_tel_description, 2: magic_tel_description}
class MAGICEventSource(EventSource): class MAGICEventSource(EventSource):
""" """
EventSource for MAGIC calibrated data. EventSource for MAGIC calibrated data.
...@@ -84,17 +97,7 @@ class MAGICEventSource(EventSource): ...@@ -84,17 +97,7 @@ class MAGICEventSource(EventSource):
# self.current_run = self._set_active_run(run_number=0) # self.current_run = self._set_active_run(run_number=0)
self.current_run = None self.current_run = None
# MAGIC telescope positions in m wrt. to the center of CTA simulations self._subarray_info = SubarrayDescription('MAGIC', magic_tel_positions, magic_tel_descriptions)
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 @staticmethod
def is_compatible(file_mask): def is_compatible(file_mask):
...@@ -229,7 +232,10 @@ class MAGICEventSource(EventSource): ...@@ -229,7 +232,10 @@ class MAGICEventSource(EventSource):
counter = 0 counter = 0
# Data container - is initialized once, and data is replaced within it after each yield # 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: # Telescopes with data:
tels_in_file = ["m1", "m2"] tels_in_file = ["m1", "m2"]
...@@ -246,6 +252,39 @@ class MAGICEventSource(EventSource): ...@@ -246,6 +252,39 @@ class MAGICEventSource(EventSource):
# Setting the new active run (class MarsRun object) # Setting the new active run (class MarsRun object)
self.current_run = self._set_active_run(run_number) 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()
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 = {1, 2}
data.mon.tel[tel_i + 1] = monitoring_camera
# Loop over the events # Loop over the events
for event_i in range(self.current_run['data'].n_stereo_events): for event_i in range(self.current_run['data'].n_stereo_events):
# Event and run ids # Event and run ids
...@@ -296,10 +335,6 @@ class MAGICEventSource(EventSource): ...@@ -296,10 +335,6 @@ class MAGICEventSource(EventSource):
# Adding event charge and peak positions per pixel # 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].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].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)
if self.is_mc == False: if self.is_mc == False:
# Adding the event arrival time # Adding the event arrival time
...@@ -345,7 +380,10 @@ class MAGICEventSource(EventSource): ...@@ -345,7 +380,10 @@ class MAGICEventSource(EventSource):
telescope = telescope.upper() telescope = telescope.upper()
# Data container - is initialized once, and data is replaced within it after each yield # 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: # Telescopes with data:
tels_in_file = ["M1", "M2"] tels_in_file = ["M1", "M2"]
...@@ -367,6 +405,38 @@ class MAGICEventSource(EventSource): ...@@ -367,6 +405,38 @@ class MAGICEventSource(EventSource):
# Setting the new active run # Setting the new active run
self.current_run = self._set_active_run(run_number) 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
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
if telescope == 'M1': if telescope == 'M1':
n_events = self.current_run['data'].n_mono_events_m1 n_events = self.current_run['data'].n_mono_events_m1
else: else:
...@@ -418,9 +488,6 @@ class MAGICEventSource(EventSource): ...@@ -418,9 +488,6 @@ class MAGICEventSource(EventSource):
# Adding event charge and peak positions per pixel # Adding event charge and peak positions per pixel
data.dl1.tel[tel_i + 1].image = event_data['image'] 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].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)
if self.is_mc == False: if self.is_mc == False:
# Adding the event arrival time # Adding the event arrival time
...@@ -465,7 +532,7 @@ class MAGICEventSource(EventSource): ...@@ -465,7 +532,7 @@ class MAGICEventSource(EventSource):
telescope = telescope.upper() telescope = telescope.upper()
# Data container - is initialized once, and data is replaced within it after each yield # Data container - is initialized once, and data is replaced within it after each yield
data = DataContainer() data = EventAndMonDataContainer()
# Telescopes with data: # Telescopes with data:
tels_in_file = ["M1", "M2"] tels_in_file = ["M1", "M2"]
...@@ -487,6 +554,35 @@ class MAGICEventSource(EventSource): ...@@ -487,6 +554,35 @@ class MAGICEventSource(EventSource):
# Setting the new active run # Setting the new active run
self.current_run = self._set_active_run(run_number) self.current_run = self._set_active_run(run_number)
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
if telescope == 'M1': if telescope == 'M1':
n_events = self.current_run['data'].n_pedestal_events_m1 n_events = self.current_run['data'].n_pedestal_events_m1
else: else:
...@@ -538,9 +634,6 @@ class MAGICEventSource(EventSource): ...@@ -538,9 +634,6 @@ class MAGICEventSource(EventSource):
# Adding event charge and peak positions per pixel # Adding event charge and peak positions per pixel
data.dl1.tel[tel_i + 1].image = event_data['image'] 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].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 # Adding the event arrival time
time_tmp = Time(event_data['mjd'], scale='utc', format='mjd') time_tmp = Time(event_data['mjd'], scale='utc', format='mjd')
...@@ -578,6 +671,8 @@ class MarsRun: ...@@ -578,6 +671,8 @@ class MarsRun:
files satisfying run_file_mask will be used. Defaults to None. files satisfying run_file_mask will be used. Defaults to None.
""" """
self.n_camera_pixels = geom.n_pixels
self.run_file_mask = run_file_mask self.run_file_mask = run_file_mask
# Preparing the lists of M1/2 data files # Preparing the lists of M1/2 data files
...@@ -609,10 +704,19 @@ class MarsRun: ...@@ -609,10 +704,19 @@ class MarsRun:
if len(run_numbers) > 1: if len(run_numbers) > 1:
raise ValueError("Run mask corresponds to more than one run: {}".format(run_numbers)) 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 = dict()
self.event_data['M1'] = self.load_events(self.m1_file_list, self.is_mc) self.event_data['M1'] = m1_data[0]
self.event_data['M2'] = self.load_events(self.m2_file_list, self.is_mc) 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 # Detecting pedestal events
self.pedestal_ids = self._find_pedestal_events() self.pedestal_ids = self._find_pedestal_events()
...@@ -621,7 +725,6 @@ class MarsRun: ...@@ -621,7 +725,6 @@ class MarsRun:
# Detecting mono events # Detecting mono events
self.mono_ids = self._find_mono_events() self.mono_ids = self._find_mono_events()
self.n_camera_pixels = 1039
@property @property
def n_events_m1(self): def n_events_m1(self):
...@@ -652,14 +755,18 @@ class MarsRun: ...@@ -652,14 +755,18 @@ class MarsRun:
return len(self.pedestal_ids['M2']) return len(self.pedestal_ids['M2'])
@staticmethod @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 Parameters
---------- ----------
file_name: str file_name: str
Name of the MAGIC calibrated file to use. Name of the MAGIC calibrated file to use.
is_mc: boolean
Specify whether Monte Carlo (True) or data (False) events are read
n_camera_pixels: int
Number of MAGIC camera pixels (not hardcoded, but specified solely via ctapipe.instrument.CameraGeometry)
Returns Returns
------- -------
...@@ -684,12 +791,24 @@ class MarsRun: ...@@ -684,12 +791,24 @@ class MarsRun:
event_data['pointing_ra'] = scipy.array([]) event_data['pointing_ra'] = scipy.array([])
event_data['pointing_dec'] = scipy.array([]) event_data['pointing_dec'] = scipy.array([])
event_data['MJD'] = scipy.array([]) event_data['MJD'] = scipy.array([])
event_data['badpixelinfo'] = []
event_data['mars_meta'] = [] event_data['mars_meta'] = []
# run-wise meta information (same for all events) # monitoring information (updated from time to time)
mars_meta = dict() monitoring_data = dict()
monitoring_data['badpixelinfo'] = []
monitoring_data['badpixelinfoMJDrange'] = []
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] event_data['file_edges'] = [0]
degrees_per_hour = 15.0 degrees_per_hour = 15.0
...@@ -732,6 +851,18 @@ class MarsRun: ...@@ -732,6 +851,18 @@ class MarsRun:
'MReportDrive.fDec' '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: # Info only applicable for MC:
mc_list = [ mc_list = [
'MMcEvt.fEnergy', 'MMcEvt.fEnergy',
...@@ -764,59 +895,90 @@ class MarsRun: ...@@ -764,59 +895,90 @@ class MarsRun:
trigger_pattern = events[b'MTriggerPattern.fPrescaled'] trigger_pattern = events[b'MTriggerPattern.fPrescaled']
stereo_event_number = events[b'MRawEvtHeader.fStereoEvtNumber'] stereo_event_number = events[b'MRawEvtHeader.fStereoEvtNumber']
# Reading meta information: # Reading run-wise meta information (same for all events in subrun):
mars_meta = dict()
mars_meta['is_simulation'] = is_mc mars_meta['is_simulation'] = is_mc
# Reading event timing information:
if not is_mc:
event_times = input_file['Events'].arrays(time_array_list)
# Computing the event arrival time
event_mjd = event_times[b'MTime.fMjd']
event_millisec = event_times[b'MTime.fTime.fMilliSec']
event_nanosec = event_times[b'MTime.fNanoSec']
event_mjd = event_mjd + (event_millisec / 1e3 + event_nanosec / 1e9) / seconds_per_day
event_data['MJD'] = scipy.concatenate((event_data['MJD'], event_mjd))
# try to read RunHeaders tree (soft fail if not present, to pass current tests)
try: try:
meta_info = input_file['RunHeaders'].arrays(metainfo_array_list) meta_info = input_file['RunHeaders'].arrays(metainfo_array_list)
mars_meta['origin'] = "MAGIC" mars_meta['origin'] = "MAGIC"
mars_meta['input_url'] = file_name mars_meta['input_url'] = file_name
mars_meta['number'] = int(meta_info[b'MRawRunHeader.fRunNumber'][0]) mars_meta['number'] = int(meta_info[b'MRawRunHeader.fRunNumber'][0])
#mars_meta['number_subrun'] = int(meta_info[b'MRawRunHeader.fSubRunIndex'][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_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 mars_meta['source_dec'] = meta_info[b'MRawRunHeader.fSourceDEC'][0] / seconds_per_hour * u.deg
is_simulation = int(meta_info[b'MRawRunHeader.fRunType'][0]) is_mc_check = int(meta_info[b'MRawRunHeader.fRunType'][0])
if is_simulation == 0: if is_mc_check == 0:
is_simulation = False is_mc_check = False
elif is_simulation == 256: elif is_mc_check == 256:
is_simulation = True is_mc_check = True
else: else:
msg = "Run type (Data or MC) of MAGIC data file not recognised." msg = "Run type (Data or MC) of MAGIC data file not recognised."
logger.error(msg) logger.error(msg)
raise ValueError(msg) raise ValueError(msg)
if is_simulation != is_mc: if is_mc_check != is_mc:
msg = "Inconsistent run type (data or MC) between file name and runheader content." msg = "Inconsistent run type (data or MC) between file name and runheader content."
logger.error(msg) logger.error(msg)
raise ValueError(msg) raise ValueError(msg)
# Reading the info only contained in real data # Reading the info only contained in real data
if is_simulation == False: if not is_mc:
badpixelinfo = input_file['RunHeaders']['MBadPixelsCam.fArray.fInfo'].array(uproot.asjagged(uproot.asdtype(np.int32))).flatten().reshape((4, 1183), order='F') 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: # now we have 4 axes:
# 0st axis: empty (?)
# 1st axis: Unsuitable pixels # 1st axis: Unsuitable pixels
# 2nd axis: Uncalibrated pixels (says why pixel is unsuitable) # 2nd axis: Uncalibrated pixels (says why pixel is unsuitable)
# 3rd axis: Bad hardware 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 # 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. unsuitable_pix_bitinfo = badpixelinfo[1][:n_camera_pixels] # take first axis
badpixelinfo = badpixelinfo[1].astype(bool) # extract unsuitable bit:
else: unsuitable_pix = np.zeros(n_camera_pixels, dtype = np.bool)
badpixelinfo = np.zeros(1183).astype(bool) for i in range(n_camera_pixels):
unsuitable_pix[i] = int('\t{0:08b}'.format(unsuitable_pix_bitinfo[i]&0xff)[-2])
monitoring_data['badpixelinfo'].append(unsuitable_pix)
# save time interval of badpixel info:
monitoring_data['badpixelinfoMJDrange'].append([event_mjd[0], event_mjd[-1]])
except KeyError: 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_mc = False
is_simulation = False
if is_simulation == False: # try to read Pedestals tree (soft fail if not present)
event_times = input_file['Events'].arrays(time_array_list) if not is_mc:
# Computing the event arrival time try:
pedestal_info = input_file['Pedestals'].arrays(pedestal_array_list)
mjd = event_times[b'MTime.fMjd']
millisec = event_times[b'MTime.fTime.fMilliSec'] pedestal_mjd = pedestal_info[b'MTimePedestals.fMjd']
nanosec = event_times[b'MTime.fNanoSec'] pedestal_millisec = pedestal_info[b'MTimePedestals.fTime.fMilliSec']
pedestal_nanosec = pedestal_info[b'MTimePedestals.fNanoSec']
mjd = mjd + (millisec / 1e3 + nanosec / 1e9) / seconds_per_day 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])
except KeyError:
logger.warning("Pedestals tree not present in file. Cleaning algorithm may fail.")
# Reading pointing information (in units of degrees): # Reading pointing information (in units of degrees):
if 'MPointingPos.' in input_file['Events']: if 'MPointingPos.' in input_file['Events']:
...@@ -859,17 +1021,17 @@ class MarsRun: ...@@ -859,17 +1021,17 @@ class MarsRun:
drive_dec_pointing_interpolator = scipy.interpolate.interp1d(drive_mjd, drive_dec, fill_value="extrapolate") drive_dec_pointing_interpolator = scipy.interpolate.interp1d(drive_mjd, drive_dec, fill_value="extrapolate")
# Interpolating the drive pointing to the event time stamps # Interpolating the drive pointing to the event time stamps
pointing_zd = drive_zd_pointing_interpolator(mjd) pointing_zd = drive_zd_pointing_interpolator(event_mjd)
pointing_az = drive_az_pointing_interpolator(mjd) pointing_az = drive_az_pointing_interpolator(event_mjd)
pointing_ra = drive_ra_pointing_interpolator(mjd) pointing_ra = drive_ra_pointing_interpolator(event_mjd)
pointing_dec = drive_dec_pointing_interpolator(mjd)