Skip to content

GitLab

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

New version of ctapipe_io_magic reading, for MC data, the info to fill MCHeaderContainer

parent 7e22a1e8
Pipeline #89331 failed with stage
in 19 minutes and 35 seconds
Hide whitespace changes
Inline Side-by-side
ctapipe_io_magic/__init__.py
View file @ d1e6159e
......@@ -148,7 +148,7 @@ class MAGICEventSource(EventSource):
for file_path in file_list:
try:
import uproot
import uproot3 as uproot
try:
with uproot.open(file_path) as input_data:
......@@ -289,7 +289,7 @@ 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:
......@@ -333,15 +333,18 @@ class MAGICEventSource(EventSource):
data.mon.tels_with_data = {1, 2}
data.mon.tel[tel_i + 1] = monitoring_camera
#fdp-start
else:
data.mcheader.num_showers = self.current_run['data'].mcheader_data['M1']['sim_nevents'] # total, including reuse
data.mcheader.energy_range_min = (self.current_run['data'].mcheader_data['M1']['sim_emin']).to(u.TeV) # GeV->TeV
data.mcheader.energy_range_max = (self.current_run['data'].mcheader_data['M1']['sim_emax']).to(u.TeV) # GeV->TeV
data.mcheader.spectral_index = self.current_run['data'].mcheader_data['M1']['sim_eslope']
data.mcheader.max_scatter_range = (self.current_run['data'].mcheader_data['M1']['sim_max_impact']).to(u.m) # cm->m
data.mcheader.max_viewcone_radius = (self.current_run['data'].mcheader_data['M1']['sim_conesemiangle']).to(u.deg) # deg->deg
#fdp-stop
data.mcheader.num_showers = self.current_run['data'].mcheader_data['M1']['sim_nevents']
data.mcheader.energy_range_min = (self.current_run['data'].mcheader_data['M1']['sim_emin']).to(u.TeV) # GeV->TeV
data.mcheader.energy_range_max = (self.current_run['data'].mcheader_data['M1']['sim_emax']).to(u.TeV) # GeV->TeV
data.mcheader.spectral_index = self.current_run['data'].mcheader_data['M1']['sim_eslope']
data.mcheader.max_scatter_range = (self.current_run['data'].mcheader_data['M1']['sim_max_impact']).to(u.m) # cm->m
data.mcheader.max_viewcone_radius = (self.current_run['data'].mcheader_data['M1']['sim_conesemiangle']).to(u.deg)# deg->deg
if data.mcheader.max_viewcone_radius != 0.:
data.mcheader.diffuse = True
else:
data.mcheader.diffuse = False
# Loop over the events
for event_i in range(self.current_run['data'].n_stereo_events):
# Event and run ids
......@@ -512,7 +515,7 @@ class MAGICEventSource(EventSource):
data.mon.tels_with_data = tels_with_data
data.mon.tel[tel_i + 1] = monitoring_camera
#fdp start
#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.energy_range_min = (self.current_run['data'].mcheader_data[telescope]['sim_emin']).to(u.TeV) # GeV->TeV
......@@ -520,13 +523,7 @@ class MAGICEventSource(EventSource):
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
#mcheader_data['sim_nevents']=int(mc_header_info[b'MMcRunHeader.fNumEvents']) #std: 5000
#mcheader_data['sim_emin']=mc_header_info[b'MMcCorsikaRunHeader.fELowLim'] #std: 10 GeV
#mcheader_data['sim_emax']=mc_header_info[b'MMcCorsikaRunHeader.fEUppLim'] #std: 30000 GeV
#mcheader_data['sim_eslope']=mc_header_info[b'MMcCorsikaRunHeader.fSlopeSpec'] #std: -1.6
#mcheader_data['sim_max_impact']=mc_header_info[b'MMcRunHeader.fImpactMax'] #std: 35000 cm
#mcheader_data['sim_conesemiangle']=mc_header_info[b'MMcRunHeader.fRandomPointingConeSemiAngle'] #std: 2.5 deg
#fdp stop
if telescope == 'M1':
n_events = self.current_run['data'].n_mono_events_m1
......@@ -831,13 +828,11 @@ class MarsRun:
self.monitoring_data['M1'] = m1_data[1]
self.monitoring_data['M2'] = m2_data[1]
#fdp start
# Getting the run-wise MC header data
if self.is_mc:
self.mcheader_data = dict()
self.mcheader_data['M1'] = m1_data[2]
self.mcheader_data['M2'] = m2_data[2]
#fdp stop
# Detecting pedestal events
self.pedestal_ids = self._find_pedestal_events()
......@@ -895,7 +890,7 @@ class MarsRun:
"""
try:
import uproot
import uproot3 as uproot
except ImportError:
msg = "The `uproot` python module is required to access the MAGIC data"
raise ImportError(msg)
......@@ -929,11 +924,9 @@ class MarsRun:
monitoring_data['PedestalFromExtractorRndm']['Mean'] = []
monitoring_data['PedestalFromExtractorRndm']['Rms'] = []
#fdp start
if is_mc:
mcheader_data = dict()
#fdp stop
#MC Header information, dictionary always created, but filled only in case of MC run
mcheader_data = dict()
event_data['file_edges'] = [0]
degrees_per_hour = 15.0
......@@ -999,7 +992,6 @@ class MarsRun:
'MMcEvt.fCoreY'
]
#fdp-start
mcheader_list = [
#'MMcRunHeader.fNumSimulatedShowers',
'MMcRunHeader.fNumEvents',
......@@ -1010,7 +1002,6 @@ class MarsRun:
#'MMcCorsikaRunHeader.fViewconeAngles',
'MMcRunHeader.fRandomPointingConeSemiAngle' # deg
]
#fdp-end
# Metadata, currently not strictly required
metainfo_array_list = [
......@@ -1021,7 +1012,9 @@ class MarsRun:
'MRawRunHeader.fSourceDEC',
'MRawRunHeader.fTelescopeNumber']
for file_name in file_list:
#print("file_name: ",file_name, is_mc, len(file_list))
input_file = uproot.open(file_name)
......@@ -1038,16 +1031,14 @@ class MarsRun:
mars_meta['is_simulation'] = is_mc
#fdp-start
if is_mc:
mc_header_info = input_file['RunHeaders'].arrays(mcheader_list)
mcheader_data['sim_nevents']=int(mc_header_info[b'MMcRunHeader.fNumEvents']) #std: 5000
mcheader_data['sim_emin']=mc_header_info[b'MMcCorsikaRunHeader.fELowLim']*u.GeV
mcheader_data['sim_emax']=mc_header_info[b'MMcCorsikaRunHeader.fEUppLim']*u.GeV
mcheader_data['sim_eslope']=mc_header_info[b'MMcCorsikaRunHeader.fSlopeSpec'] #std: -1.6
mcheader_data['sim_max_impact']=mc_header_info[b'MMcRunHeader.fImpactMax']*u.cm
mcheader_data['sim_conesemiangle']=mc_header_info[b'MMcRunHeader.fRandomPointingConeSemiAngle']*u.deg #std: 2.5 deg, also corsika viewcone is defined by "half of the cone angle".
#fdp-end
mcheader_data['sim_nevents']=int(mc_header_info[b'MMcRunHeader.fNumEvents'][0]) #std: 5000
mcheader_data['sim_emin']=mc_header_info[b'MMcCorsikaRunHeader.fELowLim'][0]*u.GeV
mcheader_data['sim_emax']=mc_header_info[b'MMcCorsikaRunHeader.fEUppLim'][0]*u.GeV
mcheader_data['sim_eslope']=mc_header_info[b'MMcCorsikaRunHeader.fSlopeSpec'][0] #std: -1.6
mcheader_data['sim_max_impact']=mc_header_info[b'MMcRunHeader.fImpactMax'][0]*u.cm
mcheader_data['sim_conesemiangle']=mc_header_info[b'MMcRunHeader.fRandomPointingConeSemiAngle'][0]*u.deg #std: 2.5 deg, also corsika viewcone is defined by "half of the cone angle".
# Reading event timing information:
if not is_mc:
......@@ -1284,11 +1275,9 @@ class MarsRun:
monitoring_data['PedestalFromExtractor'][quantity])
monitoring_data['PedestalFromExtractorRndm'][quantity] = np.array(
monitoring_data['PedestalFromExtractorRndm'][quantity])
#fdp: added mcheader_data
if is_mc:
return event_data, monitoring_data, mcheader_data
else:
return event_data, monitoring_data
return event_data, monitoring_data, mcheader_data
def _find_pedestal_events(self):
"""
......
setup.py
View file @ d1e6159e
......@@ -10,20 +10,20 @@ with open(path.join(this_directory, 'README.md'), encoding='utf-8') as f:
setup(
name='ctapipe_io_magic',
packages=find_packages(),
version='0.2.1',
version='0.2.2',
description='ctapipe plugin for reading of the calibrated MAGIC files',
long_description=long_description,
long_description_content_type='text/markdown',
install_requires=[
'ctapipe',
'astropy',
'uproot',
'uproot3',
'numpy',
'scipy'
],
tests_require=['pytest'],
setup_requires=['pytest_runner'],
author='Ievgen Vovk',
author='Ievgen Vovk et al.',
author_email='Ievgen.Vovk@mpp.mpg.de',
license='MIT',
)
Markdown is supported
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