Coverage for hiperta_stream/scripts/reorganizer_dl1hiperta300_stream_to_dl1lstchain063.py: 17%
167 statements
« prev ^ index » next coverage.py v7.4.3, created at 2024-07-16 10:16 +0000
1# !/usr/bin/env python3
3# Everything is hardcoded to LST telescopes - because this is just a temporal work that will no longer be used
4# in the moment lstchain upgrades to V0.8
6import os
7import tables
8import argparse
9import numpy as np
10from astropy.table import Table, vstack, join
11from astropy.io.misc.hdf5 import write_table_hdf5
12import copy
13import pandas as pd
14import astropy.units as u
15from lstchain.io.io import dl1_params_lstcam_key, dl1_images_lstcam_key, add_column_table
16from lstchain.reco.disp import disp
17from lstchain.reco.utils import sky_to_camera
18import time
20def add_disp_and_mc_type_to_parameters_table(dl1_file, table_path):
21 """
22 HARDCODED function obtained from `lstchain.reco.dl0_to_dl1` because `mc_alt_tel` and `mc_az_tel` are zipped within
23 `run_array_direction`.
24 1. Reconstruct the disp parameters and source position from a DL1 parameters table and write the result in the file.
25 2. Computes mc_type from the name of the file.
28 Parameters
29 ----------
30 dl1_file: HDF5 DL1 file containing the required field in `table_path`:
31 - mc_alt
32 - mc_az
33 - mc_alt_tel
34 - mc_az_tel
36 table_path: path to the parameters table in the file
38 Returns
39 -------
40 None
41 """
42 with tables.open_file(dl1_file) as hfile:
43 run_array_dir = copy.copy(hfile.root.simulation.run_config.col('run_array_direction')[0])
44 # Remember that /telescope has been moved previously
45 focal = copy.copy(hfile.root.instrument.telescope.optics.col('equivalent_focal_length')[0])
47 df = pd.read_hdf(dl1_file, key=table_path)
50 source_pos_in_camera = sky_to_camera(df.mc_alt.values * u.rad,
51 df.mc_az.values * u.rad,
52 focal * u.m,
53 run_array_dir[1] * u.rad,
54 run_array_dir[0] * u.rad,
55 )
57 disp_parameters = disp(df.x.values * u.m,
58 df.y.values * u.m,
59 source_pos_in_camera.x,
60 source_pos_in_camera.y)
62 with tables.open_file(dl1_file, mode="a") as file:
63 tab = file.root[table_path]
64 add_column_table(tab, tables.Float32Col, 'disp_dx', disp_parameters[0].value)
65 tab = file.root[table_path]
66 add_column_table(tab, tables.Float32Col, 'disp_dy', disp_parameters[1].value)
67 tab = file.root[table_path]
68 add_column_table(tab, tables.Float32Col, 'disp_norm', disp_parameters[2].value)
69 tab = file.root[table_path]
70 add_column_table(tab, tables.Float32Col, 'disp_angle', disp_parameters[3].value)
71 tab = file.root[table_path]
72 add_column_table(tab, tables.Float32Col, 'disp_sign', disp_parameters[4])
73 tab = file.root[table_path]
74 add_column_table(tab, tables.Float32Col, 'src_x', source_pos_in_camera.x.value)
75 tab = file.root[table_path]
76 add_column_table(tab, tables.Float32Col, 'src_y', source_pos_in_camera.y.value)
77 tab = file.root[table_path]
78 add_column_table(tab, tables.Float32Col, 'mc_alt_tel', np.ones(len(df)) * run_array_dir[1])
79 tab = file.root[table_path]
80 add_column_table(tab, tables.Float32Col, 'mc_az_tel', np.ones(len(df)) * run_array_dir[0])
81 if 'gamma' in dl1_file:
82 tab = file.root[table_path]
83 add_column_table(tab, tables.Float32Col, 'mc_type', np.zeros(len(df)))
84 if 'electron' in dl1_file:
85 tab = file.root[table_path]
86 add_column_table(tab, tables.Float32Col, 'mc_type', np.ones(len(df)))
87 if 'proton' in dl1_file:
88 tab = file.root[table_path]
89 add_column_table(tab, tables.Float32Col, 'mc_type', 101*np.ones(len(df)))
92def stack_and_write_images_table(input_filename, hfile_out, node_dl1_event):
93 """
94 Stack all the `tel_00X` image tables (in case they exit) and write in the v0.6 file
96 Parameters
97 input_filename : [str] input hfile name
98 hfile_out : output File pointer
99 node_dl1_event : Output hfile (V0.6) dl1.event node pointer
100 """
101 telescope_node = node_dl1_event.telescope
103 imag_per_tels = [Table(table_img.read()) for table_img in telescope_node.image]
104 image_table = vstack(imag_per_tels)
106 for tab in telescope_node.image:
107 hfile_out.remove_node(tab)
109 # Todo change names of column `image_mask` to `` ??
111 dump_plus_copy_node_to_create_new_table(input_filename,
112 hfile_out,
113 image_table,
114 hfile_out.root.dl1.event.telescope.image,
115 newname_pointer='LST_LSTCam',
116 tmp_name='imgsTable')
119def stack_and_write_parameters_table(input_filename, hfile_out, node_dl1_event, output_mc_table_pointer):
120 """
121 Stack all the `tel_00X` parameters tables (of v0.8), change names of the columns and write the table in the
122 V0.6 (lstchain like) format
124 Parameters
125 hfile_out : output File pointer
126 node_dl1_event : Output hfile (V0.6) dl1.event node pointer
127 output_mc_table_pointer : output subarray node pointer
128 """
129 telescope_node = node_dl1_event.telescope
131 param_per_tels = [Table(table_param.read()) for table_param in telescope_node.parameters]
132 parameter_table = vstack(param_per_tels)
134 for tab in telescope_node.parameters:
135 hfile_out.remove_node(tab)
137 parameter_table.rename_column('hillas_intensity', 'intensity')
138 parameter_table.rename_column('hillas_x', 'x')
139 parameter_table.rename_column('hillas_y', 'y')
140 parameter_table.rename_column('hillas_r', 'r')
141 parameter_table.rename_column('hillas_phi', 'phi')
142 parameter_table.rename_column('hillas_length', 'length')
143 parameter_table.rename_column('hillas_width', 'width')
144 parameter_table.rename_column('hillas_psi', 'psi')
145 parameter_table.rename_column('hillas_skewness', 'skewness')
146 parameter_table.rename_column('hillas_kurtosis', 'kurtosis')
147 parameter_table.rename_column('timing_slope', 'time_gradient')
148 parameter_table.rename_column('timing_intercept', 'intercept')
149 parameter_table.rename_column('morphology_num_pixels', 'n_pixels')
150 parameter_table.rename_column('morphology_num_islands', 'n_islands')
151 parameter_table.add_column(np.log10(parameter_table['intensity']), name='log_intensity')
152 parameter_table.add_column(parameter_table['width'] / parameter_table['length'], name='wl')
154 # Param table is indeed huge - it contains all the mc_events parameters (from v0.6 !!) too
155 try:
156 mc_shower_pointer = output_mc_table_pointer.mc_shower
157 except:
158 mc_shower_pointer = None
159 if mc_shower_pointer is not None:
160 mc_event_table = Table(mc_shower_pointer.read())
161 mc_event_table.remove_column('obs_id')
162 parameter_table = join(parameter_table, mc_event_table, keys='event_id')
163 parameter_table.add_column(np.log10(parameter_table['mc_energy']), name='log_mc_energy')
165 dump_plus_copy_node_to_create_new_table(input_filename,
166 hfile_out,
167 parameter_table,
168 hfile_out.root.dl1.event.telescope.parameters,
169 newname_pointer='LST_LSTCam',
170 tmp_name='paramsTable')
173def dump_plus_copy_node_to_create_new_table(input_filename, hfile_out, astropy_table_to_copy, newparent_pointer,
174 newname_pointer, tmp_name, overwrite=False):
175 """
176 General function to write an astropy table to a temporal file, and immediately after copy it to the
177 output v0.6 hfile.
179 Parameters
180 input_filename : [ste] input hfile name
181 hfile_out : output File pointer
182 astropy_table_to_copy : astropy table to be copied
183 newparent_pointer : newparent copy_node parameter
184 newname_pointer : newname copy_node parameter
185 tmp_name : [str] flag to identify the temportal table and make it unique (necessary when simultaneous reorganizers
186 are run in the same dir)
187 overwrite : overwrite parameter of the copy_node method
188 """
189 input_filename = input_filename.split('___')[0]
190 if tmp_name == '':
191 flag_name = 'UNKNOWN'
192 else:
193 flag_name = tmp_name
195 temp_table_name = f'{input_filename}_tmp_table_reoganizer_{flag_name}.tmp_hdf5_file'
196 write_table_hdf5(astropy_table_to_copy, temp_table_name, path='/root')
197 temp_table = tables.open_file(temp_table_name, 'r')
198 hfile_out.copy_node(temp_table.root.root, newparent=newparent_pointer, newname=newname_pointer, overwrite=overwrite)
199 temp_table.close()
200 os.remove(temp_table_name)
203def rename_mc_shower_colnames(input_filename, hfile_out, event_node, output_mc_table_pointer):
204 """
205 Rename column names of the `mc_shower` table and dump the table to the v0.6 output hfile.
207 Parameters
208 input_filename : [ste] input hfile name
209 hfile_out : output File pointer
210 event_node : root.dl1.event node (of output hfile, so V0.6)
211 output_mc_table_pointer : output subarray node pointer
212 """
213 mc_shower_table = Table(event_node.subarray.mc_shower.read())
214 mc_shower_table.rename_column('true_energy', 'mc_energy')
215 mc_shower_table.rename_column('true_alt', 'mc_alt')
216 mc_shower_table.rename_column('true_az', 'mc_az')
217 mc_shower_table.rename_column('true_core_x', 'mc_core_x')
218 mc_shower_table.rename_column('true_core_y', 'mc_core_y')
219 mc_shower_table.rename_column('true_h_first_int', 'mc_h_first_int')
220 mc_shower_table.rename_column('true_x_max', 'mc_x_max')
221 mc_shower_table.rename_column('true_shower_primary_id', 'mc_shower_primary_id')
223 dump_plus_copy_node_to_create_new_table(input_filename,
224 hfile_out,
225 mc_shower_table,
226 output_mc_table_pointer,
227 newname_pointer='mc_shower',
228 tmp_name='mcShowerTable',
229 overwrite=True)
232def create_hfile_out(input_filename, outfile_name, sim_pointer08, config_pointer08, dl1_pointer08, filter_pointer08):
233 """
234 Create output hfile (lstchainv0.6 like hdf5 file)
236 Parameters
237 input_filename : [ste] input hfile name
238 outfile_name : [str] output hfile name
239 sim_pointer08 : dl1-file_v0.8_simulation pointer
240 config_pointer08 : dl1-file_v.08_configuration pointer
241 dl1_pointer08 : dl1-file_v0.8_dl1 pointer
242 filter_pointer08 : dl1-file_v0.8 filters pointer
243 """
244 hfile_out = tables.open_file(outfile_name, 'w')
245 hfile_out.create_group('/', 'simulation')
246 hfile_out.create_group('/', 'dl1')
248 if sim_pointer08 is not None:
249 # Simulation node V0.6
250 # /simulation (Group) 'Simulation information of the run'
251 # children := ['mc_event' (Table), 'run_config' (Table), 'thrown_event_distribution' (Table)]
252 hfile_out.copy_node(sim_pointer08.service.shower_distribution,
253 newparent=hfile_out.root.simulation,
254 newname='thrown_event_distribution',
255 recursive=True,
256 filters=filter_pointer08)
257 hfile_out.copy_node(config_pointer08.simulation.run,
258 newparent=hfile_out.root.simulation,
259 newname='run_config',
260 recursive=True,
261 filters=filter_pointer08)
263 # Instrument node V0.6
264 # --instrument (Group)
265 # +--telescope (Group)
266 # | +--camera (Group)
267 # +--readout_LSTCam --> copied free, it can be erase.
268 # +--geometry_LSTCAM --> To be renamed to LSTCam
269 # | `--optics (Table)
270 # `--subarray (Group)
271 # `--layout (Table)
272 instrument_node = hfile_out.copy_node(config_pointer08.instrument,
273 newparent=hfile_out.root,
274 recursive=True,
275 filters=filter_pointer08)
276 hfile_out.rename_node(instrument_node.telescope.camera.geometry_LSTCam, newname='LSTCam')
278 # dl1 node V0.6
279 # +--dl1 (Group)
280 # `--event (Group)
281 # +--telescope (Group)
282 # +--image (Group)
283 # `--parameters (Group)
284 # `--subarray (Group)
285 # +--mc_shower (Table)
286 # `--trigger (Table)
287 dl1_event_node06 = hfile_out.copy_node(dl1_pointer08.event,
288 newparent=hfile_out.root.dl1,
289 recursive=True,
290 filters=filter_pointer08)
291 # This will only happen on ctapipe, not RTA
292 # hfile_out.remove_node(dl1_event_node06.telescope.trigger) # Table stored twice, remove to avoid problems.
293 try:
294 hfile_out.rename_node(dl1_event_node06.telescope.images, newname='image')
295 except tables.exceptions.NoSuchNodeError as e:
296 print(f' ** {e} Exception: No images group in the file in dl1/event/telescope/. \n'
297 f' Check the merging configuration, indeed.')
299 try:
300 subarray_pointer = hfile_out.root.dl1.event.subarray
301 # root.dl1.event.subarray.trigger
302 except tables.exceptions.NoSuchNodeError as e:
303 print(f' ** {e} Exception: No subarray pointer in dl1/event.')
304 subarray_pointer = None
306 if sim_pointer08 is not None:
307 hfile_out.copy_node(sim_pointer08.event.subarray.shower,
308 newparent=subarray_pointer,
309 newname="mc_shower",
310 recursive=True,
311 filters=filter_pointer08)
312 if subarray_pointer is not None and sim_pointer08 is not None:
313 rename_mc_shower_colnames(input_filename,
314 hfile_out,
315 dl1_event_node06,
316 subarray_pointer
317 )
319 stack_and_write_parameters_table(input_filename,
320 hfile_out,
321 dl1_event_node06,
322 subarray_pointer
323 )
325 if 'image' in dl1_event_node06.telescope:
326 stack_and_write_images_table(input_filename,
327 hfile_out,
328 dl1_event_node06
329 )
331 hfile_out.close()
334def reorganizer_dl1hiperta_stream_to_dl1lstchain063(infile, outfile, nb_tries=100, timeout=0.1):
335 """Main function to run the reorganiser"""
337 for i in range(nb_tries):
338 try:
339 hfile = tables.open_file(infile, 'r')
340 break
341 except Exception as e:
342 time.sleep(timeout)
343 print("retry "+str(i))
344 else:
345 # try 1 last time, let tables error raise if failing
346 hfile = tables.open_file(infile, 'r')
348 # dl1 v0.8 Pointers
349 try:
350 simulation_v08 = hfile.root.simulation
351 except tables.exceptions.NoSuchNodeError as e:
352 print(f' ** {e} exception: No Simulation group found in the root group of the file.')
353 simulation_v08 = None
355 configuration_v08 = hfile.root.configuration
356 dl1_v08 = hfile.root.dl1
357 filter_v08 = hfile.filters
359 create_hfile_out(infile, outfile, simulation_v08, configuration_v08, dl1_v08, filter_v08)
361 # Add disp_* and mc_type to the parameters table.
362 if simulation_v08 is not None:
363 add_disp_and_mc_type_to_parameters_table(outfile, 'dl1/event/telescope/parameters/LST_LSTCam')
365 hfile.close()
368def main():
369 """
370 Conversion from dl1 data model (ctapipe and hiper(CTA)RTA) data model, and convert it to lstchain_v0.6 data mode.
371 """
372 parser = argparse.ArgumentParser(description="Re-organize the dl1 `standard` output file from either the "
373 "hiptecta_r1_to_dl1 or hiperta_r1_dl1 to the lstchain DL1 structure")
375 parser.add_argument('--infile', '-i',
376 type=str,
377 dest='infile',
378 help='dl1 output file of `hiperta_r0_dl1` to be converted to dl1lstchain_v060',
379 default=None
380 )
382 parser.add_argument('--outfile', '-o',
383 type=str,
384 dest='outfile',
385 help='Output filename. dl1_reorganized.h5 by default.',
386 default='./dl1v0.6_reorganized.h5'
387 )
388 parser.add_argument('--nb-file-open-tries', '-n',
389 type=int,
390 dest='nb_file_tries',
391 help="Number of time the script should try to open a file (which can fail due to filesystem)",
392 default=100
393 )
394 parser.add_argument('--open-file-timeout', '-t',
395 type=float,
396 dest='file_open_timeout',
397 help="Time to wait between two attempts opening a file, in second",
398 default=0.1
399 )
401 args = parser.parse_args()
403 reorganizer_dl1hiperta_stream_to_dl1lstchain063(args.infile, args.outfile, args.nb_file_tries, args.file_open_timeout)
406if __name__ == '__main__':
407 main()