Commit fd340ded authored by Federico Di Pierro's avatar Federico Di Pierro
Browse files

new dev-dipierro-ctapipe0.12

parent 8ec7a26e
Pipeline #116246 failed with stage
in 7 seconds
......@@ -18,17 +18,24 @@ import scipy.interpolate
from astropy import units as u
from astropy.time import Time
#FDP12: eventsource is different in 0.12 w.r.t. 0.8
from ctapipe.io.eventsource import EventSource
from ctapipe.io.datalevels import DataLevel
from ctapipe.core import Container
from ctapipe.core import Field
from ctapipe.containers import DataContainer
from ctapipe.containers import EventAndMonDataContainer
#FDP12: these containers have been redefined from 0.8 to 0.12
#from ctapipe.containers import DataContainer
#from ctapipe.containers import EventAndMonDataContainer
#FDP12: these containers have not changed from 0.8 to 0.12
from ctapipe.containers import PointingContainer, TelescopePointingContainer
from ctapipe.containers import MonitoringCameraContainer
from ctapipe.containers import PedestalContainer
#FDP12: this is the new global container
from ctapipe.containers import ArrayEventContainer
#FDP12: this is the container for the mc header info
from ctapipe.containers import SimulationConfigContainer
from ctapipe.instrument import TelescopeDescription
from ctapipe.instrument import SubarrayDescription
......@@ -52,6 +59,9 @@ MAGIC_TEL_POSITIONS = {
1: [31.80, -28.10, 0.00] * u.m,
2: [-31.80, 28.10, 0.00] * u.m
}
#FDP: it is risky if tel positions are hard-coded, it would be better to read them directly from data.
#FDP: Anyway for Real MAGIC data, the coordinates are in MAGIC frame, for MC produced by sim_tel the coordinates
#FDP: are different. Could this be a problem?
# MAGIC telescope description
OPTICS = OpticsDescription.from_name('MAGIC')
......@@ -135,7 +145,7 @@ class MAGICEventSource(EventSource):
# Checking if runt type (data/MC) is consistent:
if len(is_mc_runs) > 1:
raise ValueError(
"Loaded files contain data and MC runs. Please load only data OR Monte Carlos.")
"Loaded files contain data and MC runs. Please load only data OR Monte Carlo.")
self.is_mc = is_mc_runs[0]
# Retrieving the data level (so far HARDCODED Sorcerer)
......@@ -280,18 +290,19 @@ class MAGICEventSource(EventSource):
"""
return self._stereo_event_generator()
return self._mono_event_generator() #FDP12 before the dafualt was stereo
"""FDP12
def _stereo_event_generator(self):
"""
Stereo event generator. Yields DataContainer instances, filled
with the read event data.
Returns
-------
#fdp12 ---triple comment
#fdp12Stereo event generator. Yields DataContainer instances, filled
#fdp12with the read event data.
"""
#fdp12Returns
#fdp12 -------
#fdp12---triple comment
counter = 0
# Data container - is initialized once, and data is replaced within it after each yield
......@@ -471,6 +482,8 @@ class MAGICEventSource(EventSource):
counter += 1
return
#FDP12
"""
def _mono_event_generator(self, telescope):
"""
......@@ -491,13 +504,19 @@ class MAGICEventSource(EventSource):
telescope = telescope.upper()
# Data container - is initialized once, and data is replaced within it after each yield
if not self.is_mc:
data = EventAndMonDataContainer()
else:
data = DataContainer()
#FDP12
#FDP12 if not self.is_mc:
#FDP12 data = EventAndMonDataContainer()
#FDP12 else:
#FDP12 data = DataContainer()
#FDP12 - start
data = ArrayEventContainer()
#FDP12 - end
# Telescopes with data:
tels_in_file = ["M1", "M2"]
tels_in_file = ["M1", "M2"] # in stereo is: tels_in_file = ["m1", "m2"]
if telescope not in tels_in_file:
raise ValueError("Specified telescope {:s} is not in the allowed list {}".format(
......@@ -556,18 +575,29 @@ class MAGICEventSource(EventSource):
monitoring_camera.pedestal = pedestal_info
monitoring_camera.pixel_status = badpixel_info
data.mon.tels_with_data = tels_with_data
#FDP12 - the new container has not mon.tels_with_data
#FDP12 data.mon.tels_with_data = tels_with_data
data.trigger.tels_with_trigger = tels_with_data
data.mon.tel[tel_i + 1] = monitoring_camera
#fdp (mono version not fully tested)
#fdp (mono version not fully tested)
else:
data.mcheader.num_showers = self.current_run['data'].mcheader_data[telescope]['sim_nevents'] # total, including reuse
data.mcheader.shower_reuse = self.current_run['data'].mcheader_data[telescope]['sim_reuse']
data.mcheader.energy_range_min = (self.current_run['data'].mcheader_data[telescope]['sim_emin']).to(u.TeV) # GeV->TeV
data.mcheader.energy_range_max = (self.current_run['data'].mcheader_data[telescope]['sim_emax']).to(u.TeV) # GeV->TeV
data.mcheader.spectral_index = self.current_run['data'].mcheader_data[telescope]['sim_eslope']
data.mcheader.max_scatter_range = (self.current_run['data'].mcheader_data[telescope]['sim_max_impact']).to(u.m) # cm->m
data.mcheader.max_viewcone_radius = (self.current_run['data'].mcheader_data[telescope]['sim_conesemiangle']).to(u.deg) # deg->deg
#FDP12 data.mcheader.num_showers = self.current_run['data'].mcheader_data[telescope]['sim_nevents'] # total, including reuse
#FDP12 data.mcheader.shower_reuse = self.current_run['data'].mcheader_data[telescope]['sim_reuse']
#FDP12 data.mcheader.energy_range_min = (self.current_run['data'].mcheader_data[telescope]['sim_emin']).to(u.TeV) # GeV->TeV
#FDP12 data.mcheader.energy_range_max = (self.current_run['data'].mcheader_data[telescope]['sim_emax']).to(u.TeV) # GeV->TeV
#FDP12 data.mcheader.spectral_index = self.current_run['data'].mcheader_data[telescope]['sim_eslope']
#FDP12 data.mcheader.max_scatter_range = (self.current_run['data'].mcheader_data[telescope]['sim_max_impact']).to(u.m) # cm->m
#FDP12 data.mcheader.max_viewcone_radius = (self.current_run['data'].mcheader_data[telescope]['sim_conesemiangle']).to(u.deg) # deg->deg
#FDP12: there is a new dedicated container for the MC header infos, which is not included in the general array event container
mcheader = SimulationConfigContainer()
mcheader.num_showers = self.current_run['data'].mcheader_data[telescope]['sim_nevents'] # total, including reuse
mcheader.shower_reuse = self.current_run['data'].mcheader_data[telescope]['sim_reuse']
mcheader.energy_range_min = (self.current_run['data'].mcheader_data[telescope]['sim_emin']).to(u.TeV) # GeV->TeV
mcheader.energy_range_max = (self.current_run['data'].mcheader_data[telescope]['sim_emax']).to(u.TeV) # GeV->TeV
mcheader.spectral_index = self.current_run['data'].mcheader_data[telescope]['sim_eslope']
mcheader.max_scatter_range = (self.current_run['data'].mcheader_data[telescope]['sim_max_impact']).to(u.m) # cm->m
mcheader.max_viewcone_radius = (self.current_run['data'].mcheader_data[telescope]['sim_conesemiangle']).to(u.deg) # deg->deg
if telescope == 'M1':
......@@ -586,7 +616,8 @@ class MAGICEventSource(EventSource):
event_data = self.current_run['data'].get_mono_event_data(
event_i, telescope=telescope)
data.meta = event_data['mars_meta']
#FDP12 data.meta commented out. It seems that the new containers have not the "meta" field, and it seems not used anyway
#FDP12 data.meta = event_data['mars_meta']
# Event counter
data.count = counter
......@@ -595,15 +626,17 @@ class MAGICEventSource(EventSource):
# Setting up the R0 container
data.r0.tel.clear()
data.r0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
#FDP12: trigger_type not present in the new containers
#FDP12 data.r0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
# Setting up the R1 container
data.r1.tel.clear()
data.r1.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
#FDP12 data.r1.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
# Setting up the DL0 container
data.dl0.tel.clear()
data.dl0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
#FDP12 data.dl0.tel[tel_i + 1].trigger_type = self.current_run['data'].event_data[telescope]['trigger_pattern'][event_order_number]
# Setting up the DL1 container
data.dl1.tel.clear()
......@@ -639,25 +672,38 @@ class MAGICEventSource(EventSource):
data.trigger.time = Time(
time_tmp, format='unix', scale='utc', precision=9)
else:
data.mc.energy = event_data['true_energy'] * u.GeV
data.mc.alt = (np.pi/2 - event_data['true_zd']) * u.rad
#FDP12
#FDP12 data.mc.energy = event_data['true_energy'] * u.GeV
data.simulation.shower.energy = event_data['true_energy'] * u.GeV
#FDP12 data.mc.alt = (np.pi/2 - event_data['true_zd']) * u.rad
data.simulation.shower.alt = (np.pi/2 - event_data['true_zd']) * u.rad
# check meaning of 7deg transformation (I.Vovk)
data.mc.az = -1 * \
(event_data['true_az'] - np.deg2rad(180 - 7)) * u.rad
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
#FDP12 data.mc.az = -1 * \
#FDP12 (event_data['true_az'] - np.deg2rad(180 - 7)) * u.rad
data.simulation.shower.az = -1 * (event_data['true_az'] - np.deg2rad(180 - 7)) * u.rad
#FDP12 data.mc.shower_primary_id = 1 - \
#FDP12 event_data['true_shower_primary_id']
data.simulation.shower.shower_primary_id = 1 - event_data['true_shower_primary_id']
#FDP12 data.mc.h_first_int = event_data['true_h_first_int'] * u.cm
data.simulation.shower.h_first_int = event_data['true_h_first_int'] * u.cm
# adding a 7deg rotation between the orientation of corsika (x axis = magnetic north) and MARS (x axis = geographical north) frames
# magnetic north is 7 deg westward w.r.t. geographical north
rot_corsika = 7 *u.deg
data.mc.core_x = (event_data['true_core_x']*np.cos(rot_corsika) - event_data['true_core_y']*np.sin(rot_corsika))* u.cm
data.mc.core_y = (event_data['true_core_x']*np.sin(rot_corsika) + event_data['true_core_y']*np.cos(rot_corsika))* u.cm
#FDP12 data.mc.core_x = (event_data['true_core_x']*np.cos(rot_corsika) - event_data['true_core_y']*np.sin(rot_corsika))* u.cm
#FDP12 data.mc.core_y = (event_data['true_core_x']*np.sin(rot_corsika) + event_data['true_core_y']*np.cos(rot_corsika))* u.cm
data.simulation.shower.core_x = (event_data['true_core_x']*np.cos(rot_corsika) - event_data['true_core_y']*np.sin(rot_corsika))* u.cm
data.simulation.shower.core_y = (event_data['true_core_x']*np.sin(rot_corsika) + event_data['true_core_y']*np.cos(rot_corsika))* u.cm
# 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
#FDP12 data.r0.tels_with_data = tels_with_data
#FDP12 data.r1.tels_with_data = tels_with_data
#FDP12 data.dl0.tels_with_data = tels_with_data
data.trigger.tels_with_trigger = tels_with_data
yield data
......
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