From 78c0079defeccb6650e76acc69b123e03de41ff1 Mon Sep 17 00:00:00 2001 From: Alessio Berti Date: Fri, 3 Jun 2022 14:22:13 +0200 Subject: [PATCH 01/10] Transform to telescope frame for hillas calculation, needed for coma aberration correction. --- magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py b/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py index 87d9d5cd9..3eee05752 100644 --- a/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py +++ b/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py @@ -137,7 +137,7 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False): logger.info(f'Is simulation: {is_simulation}') subarray = event_source.subarray - camera_geom = subarray.tel[tel_id].camera.geometry + camera_geom = subarray.tel[tel_id].camera.geometry.transform_to(TelescopeFrame) if is_simulation: @@ -266,8 +266,7 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False): tel_frame = TelescopeFrame(telescope_pointing=tel_pointing) - event_coord = SkyCoord(hillas_params.x, hillas_params.y, frame=camera_frame) - event_coord = event_coord.transform_to(tel_frame) + event_coord = SkyCoord(hillas_params.fov_lon, hillas_params.fov_lat, frame=tel_frame) true_disp = angular_separation( lon1=event_coord.altaz.az, From 9223486e0ca22c9c24c6fe8c7ed0149d66da4423 Mon Sep 17 00:00:00 2001 From: Alessio Berti Date: Fri, 3 Jun 2022 17:35:08 +0200 Subject: [PATCH 02/10] Use HillasParametersContainer. --- .../scripts/lst1_magic/lst1_magic_stereo_reco.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py b/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py index 6a7f25e15..4407b73e5 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py @@ -36,7 +36,7 @@ from ctapipe.containers import ( ArrayEventContainer, ImageParametersContainer, - CameraHillasParametersContainer, + HillasParametersContainer, ) from ctapipe.instrument import SubarrayDescription from magicctapipe.utils import ( @@ -213,14 +213,14 @@ def stereo_reconstruction(input_file, output_dir, config, magic_only=False): event.pointing.tel[tel_id].altitude = u.Quantity(df_tel['pointing_alt'].iloc[0], u.rad) event.pointing.tel[tel_id].azimuth = u.Quantity(df_tel['pointing_az'].iloc[0], u.rad) - hillas_params = CameraHillasParametersContainer( + hillas_params = HillasParametersContainer( intensity=float(df_tel['intensity'].iloc[0]), - x=u.Quantity(df_tel['x'].iloc[0], u.m), - y=u.Quantity(df_tel['y'].iloc[0], u.m), - r=u.Quantity(df_tel['r'].iloc[0], u.m), + fov_lon=u.Quantity(df_tel['fov_lon'].iloc[0], u.deg), + fov_lat=u.Quantity(df_tel['fov_lat'].iloc[0], u.deg), + r=u.Quantity(df_tel['r'].iloc[0], u.deg), phi=Angle(df_tel['phi'].iloc[0], u.deg), - length=u.Quantity(df_tel['length'].iloc[0], u.m), - width=u.Quantity(df_tel['width'].iloc[0], u.m), + length=u.Quantity(df_tel['length'].iloc[0], u.deg), + width=u.Quantity(df_tel['width'].iloc[0], u.deg), psi=Angle(df_tel['psi'].iloc[0], u.deg), skewness=float(df_tel['skewness'].iloc[0]), kurtosis=float(df_tel['kurtosis'].iloc[0]), From a4ad8dad3ac355bde8037beac7e152fca9239c0c Mon Sep 17 00:00:00 2001 From: Alessio Berti Date: Mon, 13 Jun 2022 10:59:02 +0200 Subject: [PATCH 03/10] Choose effective focal length for sim-telarray MC. --- magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py index 28dc10e36..697bf7442 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py @@ -250,7 +250,11 @@ def mc_dl0_to_dl1(input_file, output_dir, config, muons_analysis): n_events_skipped = 0 n_events_processed = 0 - event_source_allowed_tels = EventSource(input_file, allowed_tels=list(mc_tel_ids.values())) + event_source_allowed_tels = EventSource( + input_file, + allowed_tels=list(mc_tel_ids.values()), + focal_length_choice="effective" + ) for event in event_source_allowed_tels: From d2d8de8f70ff201900e8c13455eb5250a76d94be Mon Sep 17 00:00:00 2001 From: Alessio Berti Date: Mon, 13 Jun 2022 11:08:53 +0200 Subject: [PATCH 04/10] Construct the right container depending on frame used. --- .../lst1_magic/lst1_magic_stereo_reco.py | 42 +++++++++++++------ 1 file changed, 30 insertions(+), 12 deletions(-) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py b/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py index 4407b73e5..a2a8436cc 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py @@ -37,6 +37,7 @@ ArrayEventContainer, ImageParametersContainer, HillasParametersContainer, + CameraHillasParametersContainer, ) from ctapipe.instrument import SubarrayDescription from magicctapipe.utils import ( @@ -213,18 +214,35 @@ def stereo_reconstruction(input_file, output_dir, config, magic_only=False): event.pointing.tel[tel_id].altitude = u.Quantity(df_tel['pointing_alt'].iloc[0], u.rad) event.pointing.tel[tel_id].azimuth = u.Quantity(df_tel['pointing_az'].iloc[0], u.rad) - hillas_params = HillasParametersContainer( - intensity=float(df_tel['intensity'].iloc[0]), - fov_lon=u.Quantity(df_tel['fov_lon'].iloc[0], u.deg), - fov_lat=u.Quantity(df_tel['fov_lat'].iloc[0], u.deg), - r=u.Quantity(df_tel['r'].iloc[0], u.deg), - phi=Angle(df_tel['phi'].iloc[0], u.deg), - length=u.Quantity(df_tel['length'].iloc[0], u.deg), - width=u.Quantity(df_tel['width'].iloc[0], u.deg), - psi=Angle(df_tel['psi'].iloc[0], u.deg), - skewness=float(df_tel['skewness'].iloc[0]), - kurtosis=float(df_tel['kurtosis'].iloc[0]), - ) + if "fov_lon" in df_tel.columns: + + hillas_params = HillasParametersContainer( + intensity=float(df_tel['intensity'].iloc[0]), + fov_lon=u.Quantity(df_tel['fov_lon'].iloc[0], u.deg), + fov_lat=u.Quantity(df_tel['fov_lat'].iloc[0], u.deg), + r=u.Quantity(df_tel['r'].iloc[0], u.deg), + phi=Angle(df_tel['phi'].iloc[0], u.deg), + length=u.Quantity(df_tel['length'].iloc[0], u.deg), + width=u.Quantity(df_tel['width'].iloc[0], u.deg), + psi=Angle(df_tel['psi'].iloc[0], u.deg), + skewness=float(df_tel['skewness'].iloc[0]), + kurtosis=float(df_tel['kurtosis'].iloc[0]), + ) + + else: + + hillas_params = CameraHillasParametersContainer( + intensity=float(df_tel['intensity'].iloc[0]), + x=u.Quantity(df_tel['x'].iloc[0], u.m), + y=u.Quantity(df_tel['y'].iloc[0], u.m), + r=u.Quantity(df_tel['r'].iloc[0], u.m), + phi=Angle(df_tel['phi'].iloc[0], u.deg), + length=u.Quantity(df_tel['length'].iloc[0], u.m), + width=u.Quantity(df_tel['width'].iloc[0], u.m), + psi=Angle(df_tel['psi'].iloc[0], u.deg), + skewness=float(df_tel['skewness'].iloc[0]), + kurtosis=float(df_tel['kurtosis'].iloc[0]), + ) event.dl1.tel[tel_id].parameters = ImageParametersContainer(hillas=hillas_params) From dbc55f68cb76b25509d4a964f50d0bf23a04bcc4 Mon Sep 17 00:00:00 2001 From: YoshikiOhtani Date: Tue, 14 Jun 2022 13:45:44 +0000 Subject: [PATCH 05/10] Use the effective focal length for calculating LST length and width --- .../lst1_magic/lst1_magic_event_coincidence.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py b/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py index fab5288da..073cf9cf6 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py @@ -57,6 +57,9 @@ accuracy_time = 1e-7 # final digit of a timestamp, unit: [sec] +nominal_foclen_lst = 28 # unit: [m] +effective_foclen_lst = 29.30565 # unit: [m] + tel_names = { 1: 'LST-1', 2: 'MAGIC-I', @@ -98,10 +101,10 @@ def load_lst_data_file(input_file): try: # Try to load DL2 data at first: - event_data = pd.read_hdf(input_file, key=f'dl2/event/telescope/parameters/LST_LSTCam') + event_data = pd.read_hdf(input_file, key='dl2/event/telescope/parameters/LST_LSTCam') except: # Load DL1 data: - event_data = pd.read_hdf(input_file, key=f'dl1/event/telescope/parameters/LST_LSTCam') + event_data = pd.read_hdf(input_file, key='dl1/event/telescope/parameters/LST_LSTCam') event_data.set_index(['obs_id', 'event_id', 'tel_id'], inplace=True) event_data.sort_index(inplace=True) @@ -147,8 +150,14 @@ def load_lst_data_file(input_file): optics = pd.read_hdf(input_file, key='configuration/instrument/telescope/optics') focal_length = optics['equivalent_focal_length'][0] - event_data['length'] = focal_length * np.tan(np.deg2rad(event_data['length'])) # [deg] -> [m] - event_data['width'] = focal_length * np.tan(np.deg2rad(event_data['width'])) # [deg] -> [m] + if focal_length == nominal_foclen_lst: + + length_meter = focal_length * np.tan(np.deg2rad(event_data['length'])) + width_meter = focal_length * np.tan(np.deg2rad(event_data['width'])) + + # Convert to degrees with the effective focal length: + event_data['length'] = np.rad2deg(np.arctan2(length_meter, effective_foclen_lst)) + event_data['width'] = np.rad2deg(np.arctan2(width_meter, effective_foclen_lst)) event_data['phi'] = np.rad2deg(event_data['phi']) # [rad] -> [deg] event_data['psi'] = np.rad2deg(event_data['psi']) # [rad] -> [deg] From 00d34cb187cd0e3e64973657aac75ea2264518d5 Mon Sep 17 00:00:00 2001 From: YoshikiOhtani Date: Tue, 14 Jun 2022 17:31:22 +0000 Subject: [PATCH 06/10] Modified the LST-1 SubarrayDescription with the effective focal length --- .../scripts/lst1_magic/lst1_magic_event_coincidence.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py b/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py index 073cf9cf6..165dd94d5 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py @@ -41,6 +41,7 @@ from astropy.time import Time from ctapipe.containers import EventType from ctapipe.instrument import SubarrayDescription +from ctapipe.coordinates import CameraFrame from lstchain.reco.utils import add_delta_t_key from magicctapipe.utils import ( check_tel_combination, @@ -146,6 +147,9 @@ def load_lst_data_file(input_file): inplace=True, ) + # Read the subarray description: + subarray = SubarrayDescription.from_hdf(input_file) + # Change the units of some parameters: optics = pd.read_hdf(input_file, key='configuration/instrument/telescope/optics') focal_length = optics['equivalent_focal_length'][0] @@ -159,12 +163,12 @@ def load_lst_data_file(input_file): event_data['length'] = np.rad2deg(np.arctan2(length_meter, effective_foclen_lst)) event_data['width'] = np.rad2deg(np.arctan2(width_meter, effective_foclen_lst)) + subarray.tel[1].optics.equivalent_focal_length = u.Quantity(effective_foclen_lst, u.m) + subarray.tel[1].camera.geometry.frame = CameraFrame(focal_length=u.Quantity(effective_foclen_lst, u.m)) + event_data['phi'] = np.rad2deg(event_data['phi']) # [rad] -> [deg] event_data['psi'] = np.rad2deg(event_data['psi']) # [rad] -> [deg] - # Read the subarray description: - subarray = SubarrayDescription.from_hdf(input_file) - return event_data, subarray From 5cfbc27ff73cb031844fd61c118e6f51b7c24f65 Mon Sep 17 00:00:00 2001 From: Alessio Berti Date: Wed, 15 Jun 2022 11:14:54 +0200 Subject: [PATCH 07/10] Choose effective focal length from EventSource used for subarray. --- magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py index 697bf7442..77f9d8902 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py @@ -124,7 +124,7 @@ def mc_dl0_to_dl1(input_file, output_dir, config, muons_analysis): logger.info('\nLoading the input file:') logger.info(input_file) - event_source = EventSource(input_file) + event_source = EventSource(input_file, focal_length_choice="effective") obs_id = event_source.obs_ids[0] subarray = event_source.subarray From 360f6c1cc21bbc15e342069c9110e4e873e3f9ba Mon Sep 17 00:00:00 2001 From: YoshikiOhtani Date: Fri, 17 Jun 2022 04:21:27 +0000 Subject: [PATCH 08/10] Computes the Hillas parameters in the CameraFrame, use the effective focal length for the stereo reconstructions --- .../lst1_magic_event_coincidence.py | 12 ++---- .../lst1_magic/lst1_magic_mc_dl0_to_dl1.py | 8 +--- .../lst1_magic/lst1_magic_stereo_reco.py | 42 ++++++------------- .../scripts/lst1_magic/magic_calib_to_dl1.py | 2 +- 4 files changed, 19 insertions(+), 45 deletions(-) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py b/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py index 165dd94d5..886f5eee0 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py @@ -154,15 +154,11 @@ def load_lst_data_file(input_file): optics = pd.read_hdf(input_file, key='configuration/instrument/telescope/optics') focal_length = optics['equivalent_focal_length'][0] - if focal_length == nominal_foclen_lst: - - length_meter = focal_length * np.tan(np.deg2rad(event_data['length'])) - width_meter = focal_length * np.tan(np.deg2rad(event_data['width'])) - - # Convert to degrees with the effective focal length: - event_data['length'] = np.rad2deg(np.arctan2(length_meter, effective_foclen_lst)) - event_data['width'] = np.rad2deg(np.arctan2(width_meter, effective_foclen_lst)) + event_data['length'] = focal_length * np.tan(np.deg2rad(event_data['length'])) # [deg] -> [m] + event_data['width'] = focal_length * np.tan(np.deg2rad(event_data['width'])) # [deg] -> [m] + if focal_length == nominal_foclen_lst: + # Set the effective focal length to the subarray: subarray.tel[1].optics.equivalent_focal_length = u.Quantity(effective_foclen_lst, u.m) subarray.tel[1].camera.geometry.frame = CameraFrame(focal_length=u.Quantity(effective_foclen_lst, u.m)) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py index 77f9d8902..fe757f62f 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py @@ -124,7 +124,7 @@ def mc_dl0_to_dl1(input_file, output_dir, config, muons_analysis): logger.info('\nLoading the input file:') logger.info(input_file) - event_source = EventSource(input_file, focal_length_choice="effective") + event_source = EventSource(input_file, focal_length_choice='effective') obs_id = event_source.obs_ids[0] subarray = event_source.subarray @@ -250,11 +250,7 @@ def mc_dl0_to_dl1(input_file, output_dir, config, muons_analysis): n_events_skipped = 0 n_events_processed = 0 - event_source_allowed_tels = EventSource( - input_file, - allowed_tels=list(mc_tel_ids.values()), - focal_length_choice="effective" - ) + event_source_allowed_tels = EventSource(input_file, allowed_tels=list(mc_tel_ids.values())) for event in event_source_allowed_tels: diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py b/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py index a2a8436cc..6a7f25e15 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py @@ -36,7 +36,6 @@ from ctapipe.containers import ( ArrayEventContainer, ImageParametersContainer, - HillasParametersContainer, CameraHillasParametersContainer, ) from ctapipe.instrument import SubarrayDescription @@ -214,35 +213,18 @@ def stereo_reconstruction(input_file, output_dir, config, magic_only=False): event.pointing.tel[tel_id].altitude = u.Quantity(df_tel['pointing_alt'].iloc[0], u.rad) event.pointing.tel[tel_id].azimuth = u.Quantity(df_tel['pointing_az'].iloc[0], u.rad) - if "fov_lon" in df_tel.columns: - - hillas_params = HillasParametersContainer( - intensity=float(df_tel['intensity'].iloc[0]), - fov_lon=u.Quantity(df_tel['fov_lon'].iloc[0], u.deg), - fov_lat=u.Quantity(df_tel['fov_lat'].iloc[0], u.deg), - r=u.Quantity(df_tel['r'].iloc[0], u.deg), - phi=Angle(df_tel['phi'].iloc[0], u.deg), - length=u.Quantity(df_tel['length'].iloc[0], u.deg), - width=u.Quantity(df_tel['width'].iloc[0], u.deg), - psi=Angle(df_tel['psi'].iloc[0], u.deg), - skewness=float(df_tel['skewness'].iloc[0]), - kurtosis=float(df_tel['kurtosis'].iloc[0]), - ) - - else: - - hillas_params = CameraHillasParametersContainer( - intensity=float(df_tel['intensity'].iloc[0]), - x=u.Quantity(df_tel['x'].iloc[0], u.m), - y=u.Quantity(df_tel['y'].iloc[0], u.m), - r=u.Quantity(df_tel['r'].iloc[0], u.m), - phi=Angle(df_tel['phi'].iloc[0], u.deg), - length=u.Quantity(df_tel['length'].iloc[0], u.m), - width=u.Quantity(df_tel['width'].iloc[0], u.m), - psi=Angle(df_tel['psi'].iloc[0], u.deg), - skewness=float(df_tel['skewness'].iloc[0]), - kurtosis=float(df_tel['kurtosis'].iloc[0]), - ) + hillas_params = CameraHillasParametersContainer( + intensity=float(df_tel['intensity'].iloc[0]), + x=u.Quantity(df_tel['x'].iloc[0], u.m), + y=u.Quantity(df_tel['y'].iloc[0], u.m), + r=u.Quantity(df_tel['r'].iloc[0], u.m), + phi=Angle(df_tel['phi'].iloc[0], u.deg), + length=u.Quantity(df_tel['length'].iloc[0], u.m), + width=u.Quantity(df_tel['width'].iloc[0], u.m), + psi=Angle(df_tel['psi'].iloc[0], u.deg), + skewness=float(df_tel['skewness'].iloc[0]), + kurtosis=float(df_tel['kurtosis'].iloc[0]), + ) event.dl1.tel[tel_id].parameters = ImageParametersContainer(hillas=hillas_params) diff --git a/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py b/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py index 8322ad14b..c3f4c85c9 100644 --- a/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py +++ b/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py @@ -138,7 +138,7 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False): logger.info(f'Is simulation: {is_simulation}') subarray = event_source.subarray - camera_geom = subarray.tel[tel_id].camera.geometry.transform_to(TelescopeFrame) + camera_geom = subarray.tel[tel_id].camera.geometry if is_simulation: From 0ddbce37adfd3a3299ea419246b20efbe0679e41 Mon Sep 17 00:00:00 2001 From: YoshikiOhtani Date: Fri, 17 Jun 2022 04:26:16 +0000 Subject: [PATCH 09/10] Fix the event coordinate definition --- magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py b/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py index c3f4c85c9..0cc912097 100644 --- a/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py +++ b/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py @@ -269,7 +269,8 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False): tel_frame = TelescopeFrame(telescope_pointing=tel_pointing) - event_coord = SkyCoord(hillas_params.fov_lon, hillas_params.fov_lat, frame=tel_frame) + event_coord = SkyCoord(hillas_params.x, hillas_params.y, frame=camera_frame) + event_coord = event_coord.transform_to(tel_frame) true_disp = angular_separation( lon1=event_coord.altaz.az, From 96137c974e36849db2c4d9dc12d206345ae9cd54 Mon Sep 17 00:00:00 2001 From: YoshikiOhtani Date: Sat, 18 Jun 2022 03:11:05 +0000 Subject: [PATCH 10/10] Refactor the code --- .../lst1_magic/lst1_magic_event_coincidence.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py b/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py index 886f5eee0..2dac67f2f 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py @@ -147,9 +147,6 @@ def load_lst_data_file(input_file): inplace=True, ) - # Read the subarray description: - subarray = SubarrayDescription.from_hdf(input_file) - # Change the units of some parameters: optics = pd.read_hdf(input_file, key='configuration/instrument/telescope/optics') focal_length = optics['equivalent_focal_length'][0] @@ -157,14 +154,17 @@ def load_lst_data_file(input_file): event_data['length'] = focal_length * np.tan(np.deg2rad(event_data['length'])) # [deg] -> [m] event_data['width'] = focal_length * np.tan(np.deg2rad(event_data['width'])) # [deg] -> [m] + event_data['phi'] = np.rad2deg(event_data['phi']) # [rad] -> [deg] + event_data['psi'] = np.rad2deg(event_data['psi']) # [rad] -> [deg] + + # Read the subarray description: + subarray = SubarrayDescription.from_hdf(input_file) + if focal_length == nominal_foclen_lst: # Set the effective focal length to the subarray: subarray.tel[1].optics.equivalent_focal_length = u.Quantity(effective_foclen_lst, u.m) subarray.tel[1].camera.geometry.frame = CameraFrame(focal_length=u.Quantity(effective_foclen_lst, u.m)) - event_data['phi'] = np.rad2deg(event_data['phi']) # [rad] -> [deg] - event_data['psi'] = np.rad2deg(event_data['psi']) # [rad] -> [deg] - return event_data, subarray