Merge branch 'dev-dipierro-adapting-to-pyIRF' into 'master'

Add general run-wise MC info

See merge request !12

README.md

@@ -85,8 +85,12 @@ Monitoring data are saved in `run['data'].monitoring_data` and can also accessed
 
 Dead pixel and pedestal information are read by `magic-cta-pipe` `MAGIC_Badpixels.py` class.
 
+##### MC Header data
+Some general information about the simulated data, useful for IRF calculation, are read from the root files and stored in data.mcheader container.
+
 #### Changelog
 
 - v0.1: Initial version
 - v0.2.0: Unification of data and MC reading
 - v0.2.1: Monitoring data (Dead pixel and pedestal information)
+- v0.2.2: Added MC Header info
\ No newline at end of file

ctapipe_io_magic/__init__.py

@@ -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:
@@ -333,6 +333,19 @@ class MAGICEventSource(EventSource):
 
                     data.mon.tels_with_data = {1, 2}
                     data.mon.tel[tel_i + 1] = monitoring_camera
+            else:
+                assert self.current_run['data'].mcheader_data['M1'] == self.current_run['data'].mcheader_data['M2'], "Simulation configurations are different for M1 and M2 !!!"
+                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):
@@ -504,6 +517,16 @@ class MAGICEventSource(EventSource):
                 data.mon.tels_with_data = tels_with_data
                 data.mon.tel[tel_i + 1] = monitoring_camera
 
+             #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
+                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
+
+
             if telescope == 'M1':
                 n_events = self.current_run['data'].n_mono_events_m1
             else:
@@ -807,6 +830,12 @@ class MarsRun:
         self.monitoring_data['M1'] = m1_data[1]
         self.monitoring_data['M2'] = m2_data[1]
 
+        # 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]
+
         # Detecting pedestal events
         self.pedestal_ids = self._find_pedestal_events()
         # Detecting stereo events
@@ -863,7 +892,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)
@@ -897,6 +926,9 @@ class MarsRun:
         monitoring_data['PedestalFromExtractorRndm']['Mean'] = []
         monitoring_data['PedestalFromExtractorRndm']['Rms'] = []
 
+        #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
@@ -959,7 +991,18 @@ class MarsRun:
             'MMcEvt.fPartId',
             'MMcEvt.fZFirstInteraction',
             'MMcEvt.fCoreX',
-            'MMcEvt.fCoreY',
+            'MMcEvt.fCoreY'
+        ]
+        
+        mcheader_list = [
+            #'MMcRunHeader.fNumSimulatedShowers',
+            'MMcRunHeader.fNumEvents',
+            'MMcCorsikaRunHeader.fELowLim', #GeV
+            'MMcCorsikaRunHeader.fEUppLim', #GeV
+            'MMcCorsikaRunHeader.fSlopeSpec',
+            'MMcRunHeader.fImpactMax', #cm
+            #'MMcCorsikaRunHeader.fViewconeAngles',
+            'MMcRunHeader.fRandomPointingConeSemiAngle' # deg
         ]
 
         # Metadata, currently not strictly required
@@ -988,6 +1031,15 @@ class MarsRun:
 
             mars_meta['is_simulation'] = is_mc
 
+            if is_mc:
+                mc_header_info = input_file['RunHeaders'].arrays(mcheader_list)
+                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:
                 event_times = input_file['Events'].arrays(time_array_list)
@@ -1224,7 +1276,8 @@ class MarsRun:
                 monitoring_data['PedestalFromExtractorRndm'][quantity] = np.array(
                     monitoring_data['PedestalFromExtractorRndm'][quantity])
         
-        return event_data, monitoring_data
+        return event_data, monitoring_data, mcheader_data
+       
 
     def _find_pedestal_events(self):
         """

setup.py

@@ -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',
 )