Coverage for rta_to_lstchain_formator/reorganizer_dl1hiperta300_stream_to_dl1lstchain063.py: 61%
177 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-25 05:21 +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 argparse
7import copy
8import os
9import time
10from pathlib import Path
12import astropy.units as u
13import numpy as np
14import pandas as pd
15import tables
16from astropy.io.misc.hdf5 import write_table_hdf5
17from astropy.table import Table, join, vstack
19from rta_to_lstchain_formator.disp import disp
20from rta_to_lstchain_formator.utils.coordinates import sky_to_camera
21from rta_to_lstchain_formator.utils.pytables import add_column_table
22from rta_to_lstchain_formator.utils.rta_argparse import EnumAction
23from rta_to_lstchain_formator.utils.rta_logging import LOGGING_LEVELS, init_logging
26def add_disp_and_mc_type_to_parameters_table(dl1_file, table_path):
27 """
28 HARDCODED function obtained from `lstchain.reco.dl0_to_dl1` because `mc_alt_tel` and `mc_az_tel` are zipped within
29 `run_array_direction`.
30 1. Reconstruct the disp parameters and source position from a DL1 parameters table and write the result in the file.
31 2. Computes mc_type from the name of the file.
34 Parameters
35 ----------
36 dl1_file: HDF5 DL1 file containing the required field in `table_path`:
37 - mc_alt
38 - mc_az
39 - mc_alt_tel
40 - mc_az_tel
42 table_path: path to the parameters table in the file
44 Returns
45 -------
46 None
47 """
48 with tables.open_file(dl1_file) as hfile:
49 run_array_dir = copy.copy(hfile.root.simulation.run_config.col("run_array_direction")[0])
50 # Remember that /telescope has been moved previously
51 focal = copy.copy(hfile.root.instrument.telescope.optics.col("equivalent_focal_length")[0])
53 df = pd.read_hdf(dl1_file, key=table_path)
55 source_pos_in_camera = sky_to_camera(
56 df.mc_alt.values * u.rad,
57 df.mc_az.values * u.rad,
58 focal * u.m,
59 run_array_dir[1] * u.rad,
60 run_array_dir[0] * u.rad,
61 )
63 disp_parameters = disp(df.x.values * u.m, df.y.values * u.m, source_pos_in_camera.x, source_pos_in_camera.y)
65 with tables.open_file(dl1_file, mode="a") as file:
66 tab = file.root[table_path]
67 add_column_table(tab, tables.Float32Col, "disp_dx", disp_parameters[0].value)
68 tab = file.root[table_path]
69 add_column_table(tab, tables.Float32Col, "disp_dy", disp_parameters[1].value)
70 tab = file.root[table_path]
71 add_column_table(tab, tables.Float32Col, "disp_norm", disp_parameters[2].value)
72 tab = file.root[table_path]
73 add_column_table(tab, tables.Float32Col, "disp_angle", disp_parameters[3].value)
74 tab = file.root[table_path]
75 add_column_table(tab, tables.Float32Col, "disp_sign", disp_parameters[4])
76 tab = file.root[table_path]
77 add_column_table(tab, tables.Float32Col, "src_x", source_pos_in_camera.x.value)
78 tab = file.root[table_path]
79 add_column_table(tab, tables.Float32Col, "src_y", source_pos_in_camera.y.value)
80 tab = file.root[table_path]
81 add_column_table(tab, tables.Float32Col, "mc_alt_tel", np.ones(len(df)) * run_array_dir[1])
82 tab = file.root[table_path]
83 add_column_table(tab, tables.Float32Col, "mc_az_tel", np.ones(len(df)) * run_array_dir[0])
84 if "gamma" in dl1_file:
85 tab = file.root[table_path]
86 add_column_table(tab, tables.Float32Col, "mc_type", np.zeros(len(df)))
87 if "electron" in dl1_file:
88 tab = file.root[table_path]
89 add_column_table(tab, tables.Float32Col, "mc_type", np.ones(len(df)))
90 if "proton" in dl1_file:
91 tab = file.root[table_path]
92 add_column_table(tab, tables.Float32Col, "mc_type", 101 * np.ones(len(df)))
95def stack_and_write_images_table(input_filename, hfile_out, node_dl1_event):
96 """
97 Stack all the `tel_00X` image tables (in case they exit) and write in the v0.6 file
99 Parameters
100 input_filename : [str] input hfile name
101 hfile_out : output File pointer
102 node_dl1_event : Output hfile (V0.6) dl1.event node pointer
103 """
104 telescope_node = node_dl1_event.telescope
106 imag_per_tels = [Table(table_img.read()) for table_img in telescope_node.image]
107 image_table = vstack(imag_per_tels)
109 for tab in telescope_node.image:
110 hfile_out.remove_node(tab)
112 # Todo change names of column `image_mask` to `` ??
114 dump_plus_copy_node_to_create_new_table(
115 input_filename,
116 hfile_out,
117 image_table,
118 hfile_out.root.dl1.event.telescope.image,
119 newname_pointer="LST_LSTCam",
120 tmp_name="imgsTable",
121 )
124def stack_and_write_parameters_table(input_filename, hfile_out, node_dl1_event, output_mc_table_pointer):
125 """
126 Stack all the `tel_00X` parameters tables (of v0.8), change names of the columns and write the table in the
127 V0.6 (lstchain like) format
129 Parameters
130 hfile_out : output File pointer
131 node_dl1_event : Output hfile (V0.6) dl1.event node pointer
132 output_mc_table_pointer : output subarray node pointer
133 """
134 telescope_node = node_dl1_event.telescope
136 param_per_tels = [Table(table_param.read()) for table_param in telescope_node.parameters]
137 parameter_table = vstack(param_per_tels)
139 for tab in telescope_node.parameters:
140 hfile_out.remove_node(tab)
142 parameter_table.rename_column("hillas_intensity", "intensity")
143 parameter_table.rename_column("hillas_x", "x")
144 parameter_table.rename_column("hillas_y", "y")
145 parameter_table.rename_column("hillas_r", "r")
146 parameter_table.rename_column("hillas_phi", "phi")
147 parameter_table.rename_column("hillas_length", "length")
148 parameter_table.rename_column("hillas_width", "width")
149 parameter_table.rename_column("hillas_psi", "psi")
150 parameter_table.rename_column("hillas_skewness", "skewness")
151 parameter_table.rename_column("hillas_kurtosis", "kurtosis")
152 parameter_table.rename_column("timing_slope", "time_gradient")
153 parameter_table.rename_column("timing_intercept", "intercept")
154 parameter_table.rename_column("morphology_num_pixels", "n_pixels")
155 parameter_table.rename_column("morphology_num_islands", "n_islands")
156 parameter_table.add_column(np.log10(parameter_table["intensity"]), name="log_intensity")
157 parameter_table.add_column(parameter_table["width"] / parameter_table["length"], name="wl")
159 # Param table is indeed huge - it contains all the mc_events parameters (from v0.6 !!) too
160 try:
161 mc_shower_pointer = output_mc_table_pointer.mc_shower
162 except:
163 mc_shower_pointer = None
164 if mc_shower_pointer is not None: 164 ↛ 165line 164 didn't jump to line 165 because the condition on line 164 was never true
165 mc_event_table = Table(mc_shower_pointer.read())
166 mc_event_table.remove_column("obs_id")
167 parameter_table = join(parameter_table, mc_event_table, keys="event_id")
168 parameter_table.add_column(np.log10(parameter_table["mc_energy"]), name="log_mc_energy")
170 dump_plus_copy_node_to_create_new_table(
171 input_filename,
172 hfile_out,
173 parameter_table,
174 hfile_out.root.dl1.event.telescope.parameters,
175 newname_pointer="LST_LSTCam",
176 tmp_name="paramsTable",
177 )
180def dump_plus_copy_node_to_create_new_table(
181 input_filename, hfile_out, astropy_table_to_copy, newparent_pointer, newname_pointer, tmp_name, overwrite=False
182):
183 """
184 General function to write an astropy table to a temporal file, and immediately after copy it to the
185 output v0.6 hfile.
187 Parameters
188 input_filename : [ste] input hfile name
189 hfile_out : output File pointer
190 astropy_table_to_copy : astropy table to be copied
191 newparent_pointer : newparent copy_node parameter
192 newname_pointer : newname copy_node parameter
193 tmp_name : [str] flag to identify the temportal table and make it unique (necessary when simultaneous reorganizers
194 are run in the same dir)
195 overwrite : overwrite parameter of the copy_node method
196 """
197 input_filename = input_filename.split("___")[0]
198 if tmp_name == "": 198 ↛ 199line 198 didn't jump to line 199 because the condition on line 198 was never true
199 flag_name = "UNKNOWN"
200 else:
201 flag_name = tmp_name
203 temp_table_name = f"{input_filename}_tmp_table_reoganizer_{flag_name}.tmp_hdf5_file"
204 write_table_hdf5(astropy_table_to_copy, temp_table_name, path="/root")
205 temp_table = tables.open_file(temp_table_name, "r")
206 hfile_out.copy_node(
207 temp_table.root.root, newparent=newparent_pointer, newname=newname_pointer, overwrite=overwrite
208 )
209 temp_table.close()
210 os.remove(temp_table_name)
213def rename_mc_shower_colnames(input_filename, hfile_out, event_node, output_mc_table_pointer):
214 """
215 Rename column names of the `mc_shower` table and dump the table to the v0.6 output hfile.
217 Parameters
218 input_filename : [ste] input hfile name
219 hfile_out : output File pointer
220 event_node : root.dl1.event node (of output hfile, so V0.6)
221 output_mc_table_pointer : output subarray node pointer
222 """
223 mc_shower_table = Table(event_node.subarray.mc_shower.read())
224 mc_shower_table.rename_column("true_energy", "mc_energy")
225 mc_shower_table.rename_column("true_alt", "mc_alt")
226 mc_shower_table.rename_column("true_az", "mc_az")
227 mc_shower_table.rename_column("true_core_x", "mc_core_x")
228 mc_shower_table.rename_column("true_core_y", "mc_core_y")
229 mc_shower_table.rename_column("true_h_first_int", "mc_h_first_int")
230 mc_shower_table.rename_column("true_x_max", "mc_x_max")
231 mc_shower_table.rename_column("true_shower_primary_id", "mc_shower_primary_id")
233 dump_plus_copy_node_to_create_new_table(
234 input_filename,
235 hfile_out,
236 mc_shower_table,
237 output_mc_table_pointer,
238 newname_pointer="mc_shower",
239 tmp_name="mcShowerTable",
240 overwrite=True,
241 )
244def create_hfile_out(input_filename, outfile_name, sim_pointer08, config_pointer08, dl1_pointer08, filter_pointer08):
245 """
246 Create output hfile (lstchainv0.6 like hdf5 file)
248 Parameters
249 input_filename : [ste] input hfile name
250 outfile_name : [str] output hfile name
251 sim_pointer08 : dl1-file_v0.8_simulation pointer
252 config_pointer08 : dl1-file_v.08_configuration pointer
253 dl1_pointer08 : dl1-file_v0.8_dl1 pointer
254 filter_pointer08 : dl1-file_v0.8 filters pointer
255 """
256 hfile_out = tables.open_file(outfile_name, "w")
257 hfile_out.create_group("/", "simulation")
258 hfile_out.create_group("/", "dl1")
260 if sim_pointer08 is not None: 260 ↛ 264line 260 didn't jump to line 264 because the condition on line 260 was never true
261 # Simulation node V0.6
262 # /simulation (Group) 'Simulation information of the run'
263 # children := ['mc_event' (Table), 'run_config' (Table), 'thrown_event_distribution' (Table)]
264 hfile_out.copy_node(
265 sim_pointer08.service.shower_distribution,
266 newparent=hfile_out.root.simulation,
267 newname="thrown_event_distribution",
268 recursive=True,
269 filters=filter_pointer08,
270 )
271 hfile_out.copy_node(
272 config_pointer08.simulation.run,
273 newparent=hfile_out.root.simulation,
274 newname="run_config",
275 recursive=True,
276 filters=filter_pointer08,
277 )
279 # Instrument node V0.6
280 # --instrument (Group)
281 # +--telescope (Group)
282 # | +--camera (Group)
283 # +--readout_LSTCam --> copied free, it can be erase.
284 # +--geometry_LSTCAM --> To be renamed to LSTCam
285 # | `--optics (Table)
286 # `--subarray (Group)
287 # `--layout (Table)
288 configuration_node = hfile_out.copy_node(
289 config_pointer08, newparent=hfile_out.root, recursive=True, filters=filter_pointer08
290 )
291 hfile_out.rename_node(configuration_node.instrument.telescope.camera.geometry_LSTCam, newname="geometry_0")
292 hfile_out.rename_node(configuration_node.instrument.telescope.camera.readout_LSTCam, newname="readout_0")
294 # dl1 node V0.6
295 # +--dl1 (Group)
296 # `--event (Group)
297 # +--telescope (Group)
298 # +--image (Group)
299 # `--parameters (Group)
300 # `--subarray (Group)
301 # +--mc_shower (Table)
302 # `--trigger (Table)
303 dl1_event_node06 = hfile_out.copy_node(
304 dl1_pointer08.event, newparent=hfile_out.root.dl1, recursive=True, filters=filter_pointer08
305 )
306 # This will only happen on ctapipe, not RTA
307 # hfile_out.remove_node(dl1_event_node06.telescope.trigger) # Table stored twice, remove to avoid problems.
308 try:
309 hfile_out.rename_node(dl1_event_node06.telescope.images, newname="image")
310 except tables.exceptions.NoSuchNodeError as e:
311 print(
312 f" ** {e} Exception: No images group in the file in dl1/event/telescope/. \n"
313 f" Check the merging configuration, indeed."
314 )
316 try:
317 subarray_pointer = hfile_out.root.dl1.event.subarray
318 # root.dl1.event.subarray.trigger
319 except tables.exceptions.NoSuchNodeError as e:
320 print(f" ** {e} Exception: No subarray pointer in dl1/event.")
321 subarray_pointer = None
323 if sim_pointer08 is not None: 323 ↛ 324line 323 didn't jump to line 324 because the condition on line 323 was never true
324 hfile_out.copy_node(
325 sim_pointer08.event.subarray.shower,
326 newparent=subarray_pointer,
327 newname="mc_shower",
328 recursive=True,
329 filters=filter_pointer08,
330 )
331 if subarray_pointer is not None and sim_pointer08 is not None: 331 ↛ 332line 331 didn't jump to line 332 because the condition on line 331 was never true
332 rename_mc_shower_colnames(input_filename, hfile_out, dl1_event_node06, subarray_pointer)
334 stack_and_write_parameters_table(input_filename, hfile_out, dl1_event_node06, subarray_pointer)
336 if "image" in dl1_event_node06.telescope: 336 ↛ 339line 336 didn't jump to line 339 because the condition on line 336 was always true
337 stack_and_write_images_table(input_filename, hfile_out, dl1_event_node06)
339 hfile_out.close()
342def reorganizer_dl1hiperta_stream_to_dl1lstchain063(infile, outfile, nb_tries=100, timeout=0.1):
343 """Main function to run the reorganiser"""
345 for i in range(nb_tries): 345 ↛ 354line 345 didn't jump to line 354 because the loop on line 345 didn't complete
346 try:
347 hfile = tables.open_file(infile, "r")
348 break
349 except Exception as e:
350 time.sleep(timeout)
351 print("retry " + str(i))
352 else:
353 # try 1 last time, let tables error raise if failing
354 hfile = tables.open_file(infile, "r")
356 # dl1 v0.8 Pointers
357 try:
358 simulation_v08 = hfile.root.simulation
359 except tables.exceptions.NoSuchNodeError as e:
360 print(f" ** {e} exception: No Simulation group found in the root group of the file.")
361 simulation_v08 = None
363 configuration_v08 = hfile.root.configuration
364 dl1_v08 = hfile.root.dl1
365 filter_v08 = hfile.filters
367 create_hfile_out(infile, outfile, simulation_v08, configuration_v08, dl1_v08, filter_v08)
369 # Add disp_* and mc_type to the parameters table.
370 if simulation_v08 is not None: 370 ↛ 371line 370 didn't jump to line 371 because the condition on line 370 was never true
371 add_disp_and_mc_type_to_parameters_table(outfile, "dl1/event/telescope/parameters/LST_LSTCam")
373 hfile.close()
376def reorganiser_dl1hiperta300_stream_to_dl1lstchain063_arg_parser() -> argparse.ArgumentParser:
377 parser = argparse.ArgumentParser(
378 description="Re-organize the dl1 `standard` output file from either the "
379 "hiptecta_r1_to_dl1 or hiperta_r1_dl1 to the lstchain DL1 structure",
380 formatter_class=argparse.ArgumentDefaultsHelpFormatter,
381 )
382 parser.add_argument(
383 "--infile",
384 "-i",
385 type=str,
386 dest="infile",
387 help="dl1 output file of `hiperta_r0_dl1` to be converted to dl1lstchain_v060",
388 default=None,
389 )
390 parser.add_argument(
391 "--log_file",
392 "-l",
393 type=Path,
394 dest="log_file_path",
395 help="Path to a file where this process will log its execution",
396 )
397 parser.add_argument(
398 "--log_level",
399 type=LOGGING_LEVELS,
400 action=EnumAction,
401 dest="logging_level",
402 help="Logging level: all messages with a level severity above or equal the specified level will be logged.",
403 default=LOGGING_LEVELS.WARNING,
404 )
405 parser.add_argument(
406 "--nb-file-open-tries",
407 "-n",
408 type=int,
409 dest="nb_file_tries",
410 help="Number of time the script should try to open a file (which can fail due to filesystem)",
411 default=100,
412 )
413 parser.add_argument(
414 "--open-file-timeout",
415 "-t",
416 type=float,
417 dest="file_open_timeout",
418 help="Time to wait between two attempts opening a file, in second",
419 default=0.1,
420 )
421 parser.add_argument(
422 "--outfile",
423 "-o",
424 type=str,
425 dest="outfile",
426 help="Output filename. dl1_reorganized.h5 by default.",
427 default="./dl1v0.6_reorganized.h5",
428 )
429 return parser
432def main():
433 """
434 Conversion from dl1 data model (ctapipe and hiper(CTA)RTA) data model, and convert it to lstchain_v0.6 data mode.
435 """
437 parser = reorganiser_dl1hiperta300_stream_to_dl1lstchain063_arg_parser()
438 args = parser.parse_args()
440 init_logging(log_header_str="RTA_reorganiser", log_filename=args.log_file_path, log_level=args.logging_level)
442 reorganizer_dl1hiperta_stream_to_dl1lstchain063(
443 args.infile, args.outfile, args.nb_file_tries, args.file_open_timeout
444 )
447if __name__ == "__main__": 447 ↛ 448line 447 didn't jump to line 448 because the condition on line 447 was never true
448 main()