Note

Go to the end to download the full example code.

IO: TableLoader and EventSources#

This hands-on was presented in October 2023 at the DPPS Meeting at CEA Paris-Saclay (J. Hackfeld).

Introduction#

ctapipe provides basically two different ways of accessing its data products:

event-wise (EventSource)
column-wise (TableLoader)

EventSource(s):

EventSources read input files and generate ctapipe.containers.ArrayEventContainer instances when iterated over.

A new EventSource should be created for each type of event file read into ctapipe, e.g. simtel files are read by the ctapipe.io.SimTelEventSource .

EventSource provides a common high-level interface for accessing event information from different data sources. Creating an EventSource for a new file format or other event source ensures that data can be accessed in a common way, regardless of the file format or data origin.

EventSource itself is an abstract class, but will create an appropriate subclass if a compatible source is found for the given input_url .

TableLoader:

Loads telescope-event or subarray-event data from ctapipe HDF5 files and returns astropy.table . See Astropy docs or this video from A.Donath at the ESCAPE Summer School 2021.

This class provides high-level access to data stored in ctapipe HDF5 files, such as created by the ctapipe-process tool ( ctapipe.tools.process.ProcessorTool ).

There are multiple TableLoader methods loading data from all relevant tables (depending on the options) and joins them into single tables:

TableLoader.read_subarray_events
TableLoader.read_telescope_events
TableLoader.read_telescope_events_by_id
TableLoader.read_telescope_events_by_type

The last one returns a dict with a table per telescope type, which is needed for e.g. DL1 image data that might have different shapes for each of the telescope types as tables do not support variable length columns.

It is recommended to use the TableLoader when loading data into Tables, because its much faster than EventSources!

Code Examples#

First import some classes/modules and get the example dataset path:

 from traitlets.config import Config

 from ctapipe import utils
 from ctapipe.calib import CameraCalibrator
 from ctapipe.io import DataWriter, EventSource, TableLoader

 simtel_path = utils.get_dataset_path("gamma_prod5.simtel.zst")

EventSource(s)#

The already implemented EventSources are:

 sources = EventSource.non_abstract_subclasses()
 maxlen = max(len(source) for source in sources)
 for source in sources.values():
     print(f"{source.__name__:>{maxlen}s} -- {source.__module__}.{source.__qualname__}")

  HDF5EventSource -- ctapipe.io.hdf5eventsource.HDF5EventSource
SimTelEventSource -- ctapipe.io.simteleventsource.SimTelEventSource

EventSource will create an appropriate subclass if a compatible source is found for the given input_url :

 source = EventSource(input_url=simtel_path, max_events=5)
 print(source)

<ctapipe.io.simteleventsource.SimTelEventSource object at 0x79f08d436830>

You can now loop over the ctapipe.containers.ArrayEventContainer generated by the source.

 for event in source:
     print(event.count)

 print(repr(event))

ctapipe.containers.ArrayEventContainer:
                       index.*: event indexing information with default None
                          r0.*: Raw Data with default None
                          r1.*: R1 Calibrated Data with default None
                         dl0.*: DL0 Data Volume Reduced Data with default None
                         dl1.*: DL1 Calibrated image with default None
                         dl2.*: DL2 reconstruction info with default None
                  simulation.*: Simulated Event Information with default None
                                with type <class
                                'ctapipe.containers.SimulatedEventContainer'>
                     trigger.*: central trigger information with default None
                         count: number of events processed with default 0
                  monitoring.*: container for event-wise monitoring data (MON)
                                with default None
                        muon.*: Container for muon analysis results with default
                                None

Every time a new loop is started through the source, it tries to restart from the first event, which might not be supported by the event source. It is encouraged to use EventSource in a context manager to ensure the correct cleanups are performed when you are finished with the source:

 with EventSource(input_url=simtel_path) as source:
     for event in source:
         print(
             f"Event Count: {event.count},"
             f"Tels with trigger: {event.trigger.tels_with_trigger}"
         )

Event Count: 0,Tels with trigger: [  9  14 104]
Event Count: 1,Tels with trigger: [ 26  27  29  47  53  59  69  73 121 122 124 148 162 166]
Event Count: 2,Tels with trigger: [ 82  92 175 179]
Event Count: 3,Tels with trigger: [ 17  26  27  29  43  47  53  57  69 112 121 122 124 125 128 140 144 148
 162 166]
Event Count: 4,Tels with trigger: [  3   4   8   9  14  15  23  24  25  29  32  36  42  46  52  56 103 104
 109 110 118 119 120 126 127 129 130 137 139 143 161]
Event Count: 5,Tels with trigger: [ 87  98 151]
Event Count: 6,Tels with trigger: [ 68  72  76  94  98 151 165 177]

You can hand in the ID’s of the telescopes to be included in the data with the allowed_tels attribute. If given, only this subset of telescopes will be present in the generated events. If None, all available telescopes are used.

 eventsource_config = Config(
     {"EventSource": {"max_events": 5, "allowed_tels": [3, 4, 9]}}
 )

 with EventSource(input_url=simtel_path, config=eventsource_config) as source:
     for event in source:
         print(
             f"Event Count: {event.count},"
             f"Tels with trigger: {event.trigger.tels_with_trigger}"
         )

Event Count: 0,Tels with trigger: [9]
Event Count: 1,Tels with trigger: [3 4 9]

If you want to calibrate your data in the event loop and write it to an .h5 file with the DataWriter :

 source = EventSource(input_url=simtel_path, max_events=50)
 calibrate = CameraCalibrator(subarray=source.subarray)

 with DataWriter(
     event_source=source,
     output_path="events.dl1.h5",
     write_dl1_parameters=False,
     overwrite=True,
     write_dl1_images=True,
 ) as write_data:
     for event in source:
         calibrate(event)
         write_data(event)

Alternatively doing it with ctapipe-process would look like this:

! ctapipe-process -i {simtel_path} -o events.dl1.h5 --overwrite --progress

TableLoader#

Create a TableLoader instance with the above created dl1 file:

 loader = TableLoader(input_url="events.dl1.h5")

Alternatively using a config file:

 tableloader_config = Config(
     {
         "TableLoader": {
             "input_url": "events.dl1.h5",
         }
     }
 )

 loader = TableLoader(config=tableloader_config)

Reading subarray-based event information:

 subarray_events = loader.read_subarray_events(
     start=None,
     stop=None,
     dl2=False,
     simulated=True,
     observation_info=False,
 )

 subarray_events

Table length=7

obs_id	event_id	time	tels_with_trigger	event_type	true_energy	true_alt	true_az	true_core_x	true_core_y	true_h_first_int	true_x_max	true_starting_grammage	true_shower_primary_id
					TeV	deg	deg	m	m	m	g / cm2	g / cm2
int32	int64	Time	bool[180]	int64	float64	float64	float64	float64	float64	float64	float64	float64	int64
1	4009	1604409937.84765	False .. False	32	0.07287972420454025	69.99999967119774	359.999982697156	-146.42298889160156	341.6435852050781	32396.98828125	159.41175842285156	0.0	0
1	5101	1604409953.115399	False .. False	32	1.7454713582992554	69.99999967119774	359.999982697156	-737.0097045898438	-266.6009521484375	20947.525390625	286.7567443847656	0.0	0
1	5103	1604409960.710867	False .. False	32	1.7454713582992554	69.99999967119774	359.999982697156	791.504150390625	914.0662231445312	20947.525390625	286.7567443847656	0.0	0
1	5104	1604409965.387411	False .. False	32	1.7454713582992554	69.99999967119774	359.999982697156	-534.9799194335938	-443.40386962890625	20947.525390625	286.7567443847656	0.0	0
1	5105	1604409973.646828	False .. False	32	1.7454713582992554	69.99999967119774	359.999982697156	-324.859130859375	398.5491027832031	20947.525390625	286.7567443847656	0.0	0
1	5108	1604409990.315764	False .. False	32	1.7454713582992554	69.99999967119774	359.999982697156	-1299.1031494140625	217.68312072753906	20947.525390625	286.7567443847656	0.0	0
1	5109	1604409992.675785	False .. False	32	1.7454713582992554	69.99999967119774	359.999982697156	-994.0153198242188	485.241943359375	20947.525390625	286.7567443847656	0.0	0

Reading subarray-based event information in chunks:

 subarray_events_chunked = loader.read_subarray_events_chunked(
     chunk_size=3,
     dl2=False,
     simulated=True,
     observation_info=False,
 )

 for chunk in subarray_events_chunked:
     print(" \n", chunk)

 Chunk(start=0, stop=3, data=<Table length=3>
obs_id event_id ... true_starting_grammage true_shower_primary_id
                ...        g / cm2
int32   int64   ...        float64                 int64
------ -------- ... ---------------------- ----------------------
     1     4009 ...                    0.0                      0
     1     5101 ...                    0.0                      0
     1     5103 ...                    0.0                      0)

 Chunk(start=3, stop=6, data=<Table length=3>
obs_id event_id ... true_starting_grammage true_shower_primary_id
                ...        g / cm2
int32   int64   ...        float64                 int64
------ -------- ... ---------------------- ----------------------
     1     5104 ...                    0.0                      0
     1     5105 ...                    0.0                      0
     1     5108 ...                    0.0                      0)

 Chunk(start=6, stop=7, data=<Table length=1>
obs_id event_id ... true_starting_grammage true_shower_primary_id
                ...        g / cm2
int32   int64   ...        float64                 int64
------ -------- ... ---------------------- ----------------------
     1     5109 ...                    0.0                      0)

Reading just LST events:

 lst_events = loader.read_telescope_events(
     telescopes=[1, 2, 3, 4],
     start=None,
     stop=None,
     dl1_images=True,
     dl1_parameters=False,
     dl1_muons=False,
     dl2=False,
     simulated=False,
     true_images=True,
     true_parameters=False,
     instrument=True,
     observation_info=False,
 )

 lst_events

Table length=2

obs_id	event_id	tel_id	time_mono	event_type_mono	n_trigger_pixels	image	peak_time	is_valid	true_image_sum	true_image	name	type	pos_x	pos_y	pos_z	camera_name	optics_name_1	tel_description	optics_name_2	size_type	reflector_shape	mirror_area	n_mirrors	n_mirror_tiles	equivalent_focal_length	effective_focal_length	telescope_pointing_azimuth	telescope_pointing_altitude	time	tels_with_trigger	event_type
													m	m	m							m2			m	m	rad	rad
int32	int64	int32	Time	int64	int64	float32[1855]	float32[1855]	bool	int32	int32[1855]	str5	str3	float32	float32	float32	str9	str5	str17	str5	str3	str20	float64	int64	int64	float64	float64	float32	float32	Time	bool[180]	int64
1	5105	3	1604409973.646828	32	-1	-1.5701737 .. -0.9145646	27.381046 .. 29.717525	True	94	0 .. 0	LST	LST	-19.396	65.2	31.0	LSTCam	LST	LST_LST_LSTCam	LST	LST	PARABOLIC	386.7332458496094	1	198	28.0	29.30565071105957	0.0	1.2217305	1604409973.646828	False .. False	32
1	5105	4	1604409973.646828	32	-1	1.0622615 .. 0.9853555	2.7207987 .. 37.355812	True	175	0 .. 0	LST	LST	-120.033	1.151	33.1	LSTCam	LST	LST_LST_LSTCam	LST	LST	PARABOLIC	386.7332458496094	1	198	28.0	29.30565071105957	0.0	1.2217305	1604409973.646828	False .. False	32

Loading telescope events by type returns a dict with the different telescope types:

 telescope_events_by_type = loader.read_telescope_events_by_type(
     telescopes=["LST_LST_LSTCam", "MST_MST_FlashCam"],
     start=None,
     stop=None,
     dl1_images=True,
     dl1_parameters=False,
     dl1_muons=False,
     dl2=False,
     simulated=False,
     true_images=True,
     true_parameters=False,
     instrument=True,
     observation_info=False,
 )

 for tel_type, table in telescope_events_by_type.items():
     print(f"Telescope Type: {tel_type} \n", table, "\n")

Telescope Type: LST_LST_LSTCam
 obs_id event_id tel_id ...        time       tels_with_trigger event_type
                       ...
------ -------- ------ ... ----------------- ----------------- ----------
   5105      3 ... 1604409973.646828    False .. False         32
   5105      4 ... 1604409973.646828    False .. False         32

Telescope Type: MST_MST_FlashCam
 obs_id event_id tel_id ...        time       tels_with_trigger event_type
                       ...
------ -------- ------ ... ----------------- ----------------- ----------
   4009      9 ...  1604409937.84765    False .. False         32
   4009     14 ...  1604409937.84765    False .. False         32
   5101     26 ... 1604409953.115399    False .. False         32
   5101     27 ... 1604409953.115399    False .. False         32
   5101     29 ... 1604409953.115399    False .. False         32
   5104     17 ... 1604409965.387411    False .. False         32
   5104     26 ... 1604409965.387411    False .. False         32
   5104     27 ... 1604409965.387411    False .. False         32
   5104     29 ... 1604409965.387411    False .. False         32
   5104    125 ... 1604409965.387411    False .. False         32
   5105      8 ... 1604409973.646828    False .. False         32
   5105      9 ... 1604409973.646828    False .. False         32
   5105     14 ... 1604409973.646828    False .. False         32
   5105     15 ... 1604409973.646828    False .. False         32
   5105     23 ... 1604409973.646828    False .. False         32
   5105     24 ... 1604409973.646828    False .. False         32
   5105     25 ... 1604409973.646828    False .. False         32
   5105     29 ... 1604409973.646828    False .. False         32
   5105    126 ... 1604409973.646828    False .. False         32
   5105    127 ... 1604409973.646828    False .. False         32

Loading telescope events by ID returns a dict with the different telescope IDs:

 telescope_events_by_id = loader.read_telescope_events_by_id(
     telescopes=[3, 14],
     start=None,
     stop=None,
     dl1_images=True,
     dl1_parameters=False,
     dl1_muons=False,
     dl2=False,
     simulated=False,
     true_images=True,
     true_parameters=False,
     instrument=True,
     observation_info=False,
 )
 for tel_id, table in telescope_events_by_id.items():
     print(f"Telescope ID: {tel_id} \n", table, "\n")

Telescope ID: 3
 obs_id event_id tel_id ...        time       tels_with_trigger event_type
                       ...
------ -------- ------ ... ----------------- ----------------- ----------
     1     5105      3 ... 1604409973.646828    False .. False         32

Telescope ID: 14
 obs_id event_id tel_id ...        time       tels_with_trigger event_type
                       ...
------ -------- ------ ... ----------------- ----------------- ----------
     1     4009     14 ...  1604409937.84765    False .. False         32
     1     5105     14 ... 1604409973.646828    False .. False         32

read_telescope_events_chunked
read_telescope_events_by_type_chunked
read_telescope_events_by_id_chunked

are also available.

Reading e.g. simulation- or observation-information:

 simulation_configuration = loader.read_simulation_configuration()
 observation_information = loader.read_observation_information()

 simulation_configuration

Table length=1

obs_id	run_number	corsika_version	simtel_version	energy_range_min	energy_range_max	prod_site_B_total	prod_site_B_declination	prod_site_B_inclination	prod_site_alt	spectral_index	shower_prog_start	shower_prog_id	detector_prog_start	detector_prog_id	n_showers	shower_reuse	max_alt	min_alt	max_az	min_az	diffuse	max_viewcone_radius	min_viewcone_radius	max_scatter_range	min_scatter_range	core_pos_mode	atmosphere	corsika_iact_options	corsika_low_E_model	corsika_high_E_model	corsika_bunchsize	corsika_wlen_min	corsika_wlen_max	corsika_low_E_detail	corsika_high_E_detail
				TeV	TeV	uT	rad	rad	m								rad	rad	rad	rad		deg	deg	m	m							nm	nm
int32	int32	int64	int64	float32	float32	float64	float64	float64	float64	float64	int64	int64	int64	int64	int64	int64	float32	float32	float32	float32	int64	float32	float32	float32	float32	int64	int64	int64	int64	int64	float64	float64	float64	int64	int64
1	1	7710	1593349643	0.003	330.0	22.585954666137695	-0.07911577820777893	-0.42764246463775635	2147.0	-2.0	1604404800	1	1604409883	1	100	10	1.2217305	1.2217305	7.987575e-08	7.987575e-08	0	0.0	0.0	2000.0	0.0	1	99	187	2	3	5.0	240.0	1000.0	0	303

Now you have astropy.table s including all the relevant data you need for your analyses.

Total running time of the script: (0 minutes 6.766 seconds)

Gallery generated by Sphinx-Gallery