diff --git a/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py b/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py index 1bc091756..87d9d5cd9 100644 --- a/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py +++ b/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py @@ -2,17 +2,16 @@ # coding: utf-8 """ -Author: Yoshiki Ohtani (ICRR, ohtani@icrr.u-tokyo.ac.jp) - -This script processes events of MAGIC calibrated data (*_Y_*.root) with the MARS-like image cleaning +This script processes the events of MAGIC calibrated data (*_Y_*.root) with the MARS-like image cleaning and computes the DL1 parameters (i.e., Hillas, timing and leakage parameters). It saves only the events that all the DL1 parameters are successfully reconstructed. The telescope IDs are reset to the following ones for the combined analysis with LST-1, whose telescope ID is 1: MAGIC-I: tel_id = 2, MAGIC-II: tel_id = 3 -It searches for all the subrun files belonging to the same observation ID and stored in the same directory as an input subrun file. -Then it reads drive reports from the files and uses the information to reconstruct the telescope pointing direction. -If the "--process-run" argument is given, it not only reads drive reports but also processes all the events of the subrun files. +When an input is real data, the script searches for all the subrun files belonging to the same observation ID +and stored in the same directory as an input subrun file. Then it reads drive reports from the files and uses the information +to reconstruct the telescope pointing direction. Thus, it is best to store all the files in the same directory. +If the "--process-run" argument is given, it not only reads drive reports but also processes all the events of the subrun files at once. Usage: $ python magic_calib_to_dl1.py @@ -22,6 +21,7 @@ (--process-run) """ +import re import time import yaml import logging @@ -30,8 +30,11 @@ import numpy as np from pathlib import Path from astropy import units as u +from astropy.coordinates import AltAz, SkyCoord +from astropy.coordinates.angle_utilities import angular_separation from ctapipe.io import HDF5TableWriter from ctapipe.core import Container, Field +from ctapipe.coordinates import CameraFrame, TelescopeFrame from ctapipe.image import ( number_of_islands, hillas_parameters, @@ -41,6 +44,7 @@ from ctapipe.instrument import SubarrayDescription from ctapipe_io_magic import MAGICEventSource from magicctapipe.image import MAGICClean +from magicctapipe.utils import calculate_impact logger = logging.getLogger(__name__) logger.addHandler(logging.StreamHandler()) @@ -51,11 +55,6 @@ sec2nsec = 1e9 -tel_positions = { - 2: u.Quantity([39.3, -62.55, -0.97], u.m), # MAGIC-I - 3: u.Quantity([-31.21, -14.57, 0.2], u.m), # MAGIC-II -} - pedestal_types = [ 'fundamental', 'from_extractor', @@ -64,6 +63,7 @@ __all__ = [ 'EventInfoContainer', + 'SimEventInfoContainer', 'magic_calib_to_dl1', ] @@ -82,6 +82,25 @@ class EventInfoContainer(Container): n_islands = Field(-1, 'Number of islands of a cleaned image') +class SimEventInfoContainer(Container): + """ Container to store simulated event information """ + + obs_id = Field(-1, 'Observation ID') + event_id = Field(-1, 'Event ID') + tel_id = Field(-1, 'Telescope ID') + pointing_alt = Field(-1, 'Telescope pointing altitude', u.rad) + pointing_az = Field(-1, 'Telescope pointing azimuth', u.rad) + true_energy = Field(-1, 'Simulated event true energy', u.TeV) + true_alt = Field(-1, 'Simulated event true altitude', u.deg) + true_az = Field(-1, 'Simulated event true azimuth', u.deg) + true_disp = Field(-1, 'Simulated event true disp', u.deg) + true_core_x = Field(-1, 'Simulated event true core x', u.m) + true_core_y = Field(-1, 'Simulated event true core y', u.m) + true_impact = Field(-1, 'Simulated event true impact', u.m) + n_pixels = Field(-1, 'Number of pixels of a cleaned image') + n_islands = Field(-1, 'Number of islands of a cleaned image') + + def magic_calib_to_dl1(input_file, output_dir, config, process_run=False): """ Processes MAGIC calibrated events and computes the DL1 parameters. @@ -95,54 +114,89 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False): config: dict Configuration for the LST-1 + MAGIC analysis process_run: bool - If True, it processes events of all the subrun files - belonging to the same observation ID as an input subrun file + If True, it processes the events of all the subrun files at once """ config_cleaning = config['MAGIC']['magic_clean'] - if config_cleaning['find_hotpixels'] == 'auto': - logger.info('\nSetting the "find_hotpixels" option to True...') - config_cleaning.update({'find_hotpixels': True}) - - pedestal_type = config_cleaning.pop('pedestal_type') - - if pedestal_type not in pedestal_types: - raise KeyError(f'Unknown pedestal type "{pedestal_type}".') - logger.info('\nConfiguration for the image cleaning:') logger.info(config_cleaning) # Load the input file: + logger.info('\nLoading the input file:') + logger.info(input_file) + event_source = MAGICEventSource(input_file, process_run=process_run) obs_id = event_source.obs_ids[0] tel_id = event_source.telescope + is_simulation = event_source.is_simulation - logger.info(f'\nProcess the following data (process_run = {process_run}):') - for root_file in event_source.file_list: - logger.info(root_file) + logger.info(f'\nObservation ID: {obs_id}') + logger.info(f'Telescope ID: {tel_id}') + logger.info(f'Is simulation: {is_simulation}') - # Configure the MAGIC image cleaning: - i_ped_type = np.where(np.array(pedestal_types) == pedestal_type)[0][0] + subarray = event_source.subarray + camera_geom = subarray.tel[tel_id].camera.geometry + + if is_simulation: + + if config_cleaning['find_hotpixels'] is not False: + logger.warning('\nHot pixels do not exist in a simulation. Setting the "find_hotpixels" option to False...') + config_cleaning.update({'find_hotpixels': False}) - camera_geom = event_source.subarray.tel[tel_id].camera.geometry + unsuitable_mask = None + + tel_position = subarray.positions[tel_id] + focal_length = subarray.tel[tel_id].optics.equivalent_focal_length + + camera_frame = CameraFrame( + rotation=camera_geom.cam_rotation, + focal_length=focal_length, + ) + + else: + pedestal_type = config_cleaning.pop('pedestal_type') + + if pedestal_type not in pedestal_types: + raise KeyError(f'Unknown pedestal type "{pedestal_type}".') + + i_ped_type = np.where(np.array(pedestal_types) == pedestal_type)[0][0] + + if config_cleaning['find_hotpixels'] == 'auto': + logger.info('\nSetting the "find_hotpixels" option to True...') + config_cleaning.update({'find_hotpixels': True}) + + if process_run: + logger.info(f'\nProcess the following data (process_run = True):') + for root_file in event_source.file_list: + logger.info(root_file) + + # Configure the MAGIC image cleaning: magic_clean = MAGICClean(camera_geom, config_cleaning) # Prepare for saving data to an output file: Path(output_dir).mkdir(exist_ok=True, parents=True) - if process_run: - output_file = f'{output_dir}/dl1_M{tel_id}.Run{obs_id:08}.h5' + if is_simulation: + regex = r'GA_M\d_(\S+)_\d_\d+_Y_*' + file_name = Path(input_file).resolve().name + + parser = re.findall(regex, file_name)[0] + output_file = f'{output_dir}/dl1_M{tel_id}_GA_{parser}.Run{obs_id}.h5' + else: - subrun_id = event_source.metadata['subrun_number'][0] - output_file = f'{output_dir}/dl1_M{tel_id}.Run{obs_id:08}.{subrun_id:03}.h5' + if process_run: + output_file = f'{output_dir}/dl1_M{tel_id}.Run{obs_id:08}.h5' + else: + subrun_id = event_source.metadata['subrun_number'][0] + output_file = f'{output_dir}/dl1_M{tel_id}.Run{obs_id:08}.{subrun_id:03}.h5' # Start processing the events: - n_events_skipped = 0 - logger.info('\nProcessing the events...') + n_events_skipped = 0 + with HDF5TableWriter(output_file, group_name='events', mode='w') as writer: for event in event_source: @@ -150,15 +204,14 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False): if event.count % 100 == 0: logger.info(f'{event.count} events') - # Get bad pixels: - dead_pixels = event.mon.tel[tel_id].pixel_status.hardware_failing_pixels[0] - badrms_pixels = event.mon.tel[tel_id].pixel_status.pedestal_failing_pixels[i_ped_type] - unsuitable_mask = np.logical_or(dead_pixels, badrms_pixels) + if not is_simulation: + dead_pixels = event.mon.tel[tel_id].pixel_status.hardware_failing_pixels[0] + badrms_pixels = event.mon.tel[tel_id].pixel_status.pedestal_failing_pixels[i_ped_type] + unsuitable_mask = np.logical_or(dead_pixels, badrms_pixels) # Apply the image cleaning: signal_pixels, image, peak_time = magic_clean.clean_image(event.dl1.tel[tel_id].image, - event.dl1.tel[tel_id].peak_time, - unsuitable_mask) + event.dl1.tel[tel_id].peak_time, unsuitable_mask) image_cleaned = image.copy() image_cleaned[~signal_pixels] = 0 @@ -203,25 +256,75 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False): n_events_skipped += 1 continue - # Set the event information to the container. - # To keep the precision of a timestamp for the event coincidence with LST-1, - # here we set the integral and fractional parts separately as "time_sec" and "time_nanosec": - timestamp = event.trigger.tel[tel_id].time.to_value(format='unix', subfmt='long') - fractional, integral = np.modf(timestamp) - - time_sec = int(np.round(integral)) - time_nanosec = int(np.round(fractional * sec2nsec)) - - event_info = EventInfoContainer( - obs_id=event.index.obs_id, - event_id=event.index.event_id, - pointing_alt=event.pointing.tel[tel_id].altitude, - pointing_az=event.pointing.tel[tel_id].azimuth, - time_sec=time_sec, - time_nanosec=time_nanosec, - n_pixels=n_pixels, - n_islands=n_islands, - ) + if is_simulation: + + # Compute the DISP parameter: + tel_pointing = AltAz( + alt=event.pointing.tel[tel_id].altitude, + az=event.pointing.tel[tel_id].azimuth, + ) + + 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) + + true_disp = angular_separation( + lon1=event_coord.altaz.az, + lat1=event_coord.altaz.alt, + lon2=event.simulation.shower.az, + lat2=event.simulation.shower.alt, + ) + + # Calculate the impact parameter: + true_impact = calculate_impact( + core_x=event.simulation.shower.core_x, + core_y=event.simulation.shower.core_y, + az=event.simulation.shower.az, + alt=event.simulation.shower.alt, + tel_pos_x=tel_position[0], + tel_pos_y=tel_position[1], + tel_pos_z=tel_position[2], + ) + + # Set the event information to the container: + event_info = SimEventInfoContainer( + obs_id=event.index.obs_id, + event_id=event.index.event_id, + pointing_alt=event.pointing.tel[tel_id].altitude, + pointing_az=event.pointing.tel[tel_id].azimuth, + true_energy=event.simulation.shower.energy, + true_alt=event.simulation.shower.alt, + true_az=event.simulation.shower.az, + true_disp=true_disp, + true_core_x=event.simulation.shower.core_x, + true_core_y=event.simulation.shower.core_y, + true_impact=true_impact, + n_pixels=n_pixels, + n_islands=n_islands, + ) + + else: + timestamp = event.trigger.tel[tel_id].time.to_value(format='unix', subfmt='long') + + # To keep the precision of a timestamp for the event coincidence with LST-1, + # here we set the integral and fractional parts separately as "time_sec" and "time_nanosec": + fractional, integral = np.modf(timestamp) + + time_sec = int(np.round(integral)) + time_nanosec = int(np.round(fractional * sec2nsec)) + + # Set the event information to the container: + event_info = EventInfoContainer( + obs_id=event.index.obs_id, + event_id=event.index.event_id, + pointing_alt=event.pointing.tel[tel_id].altitude, + pointing_az=event.pointing.tel[tel_id].azimuth, + time_sec=time_sec, + time_nanosec=time_nanosec, + n_pixels=n_pixels, + n_islands=n_islands, + ) # Reset the telescope IDs: if tel_id == 1: @@ -238,17 +341,24 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False): logger.info(f'\nIn total {n_events_processed} events are processed.') logger.info(f'({n_events_skipped} events are skipped)') - # Reset the telescope IDs of the telescope descriptions: + # Reset the telescope IDs of the subarray descriptions: + tel_positions = { + 2: subarray.positions[1], # MAGIC-I + 3: subarray.positions[2], # MAGIC-II + } + tel_descriptions = { 2: event_source.subarray.tel[1], # MAGIC-I 3: event_source.subarray.tel[2], # MAGIC-II } - # Save the subarray description. - # Here we save the MAGIC telescope positions relative to the center of the LST-1 + MAGIC array, - # which are also used for sim_telarray simulations: - subarray = SubarrayDescription('MAGIC-Array', tel_positions, tel_descriptions) - subarray.to_hdf(output_file) + subarray_magic = SubarrayDescription('MAGIC-Array', tel_positions, tel_descriptions) + subarray_magic.to_hdf(output_file) + + if is_simulation: + # Save the simulation configuration: + with HDF5TableWriter(output_file, group_name='simulation', mode='a') as writer: + writer.write('config', event_source.simulation_config[obs_id]) logger.info('\nOutput file:') logger.info(output_file) @@ -277,7 +387,7 @@ def main(): parser.add_argument( '--process-run', dest='process_run', action='store_true', - help='Process all the subrun files of the same observation ID at once.', + help='Process the events of all the subrun files at once.', ) args = parser.parse_args() diff --git a/magicctapipe/scripts/lst1_magic/merge_hdf_files.py b/magicctapipe/scripts/lst1_magic/merge_hdf_files.py index 23cf93bbe..8cf25e110 100644 --- a/magicctapipe/scripts/lst1_magic/merge_hdf_files.py +++ b/magicctapipe/scripts/lst1_magic/merge_hdf_files.py @@ -139,14 +139,20 @@ def merge_hdf_files(input_dir, run_wise=False, subrun_wise=False): file_names_unique = np.unique(file_names) run_ids_unique = np.unique(run_ids) - if len(file_names_unique) == 1: - file_name = file_names_unique[0] + n_file_names = len(file_names_unique) - elif file_names_unique.tolist() == ['dl1_M1.Run', 'dl1_M2.Run']: - file_name = 'dl1_MAGIC.Run' + if n_file_names == 1: + file_name = file_names_unique[0] else: - raise RuntimeError('Multiple types of files are found in the input directory.') + replaced_name = file_names_unique[0].replace('M1', 'M2') + is_same_type = (file_names_unique[1] == replaced_name) + + if (n_file_names == 2) and is_same_type: + file_name = file_names_unique[0].replace('M1', 'MAGIC') + + else: + RuntimeError('Multiple types of files are found in the input directory.') # Merge the input files: output_dir = f'{input_dir}/merged' diff --git a/magicctapipe/utils/functions.py b/magicctapipe/utils/functions.py index 6079be7c2..ed76e09dc 100644 --- a/magicctapipe/utils/functions.py +++ b/magicctapipe/utils/functions.py @@ -215,8 +215,8 @@ def transform_altaz_to_radec(alt, az, timestamp): """ # Hardcode the longitude/latitude of ORM: - LAT_ORM = u.Quantity(28.76177, u.deg) LON_ORM = u.Quantity(-17.89064, u.deg) + LAT_ORM = u.Quantity(28.76177, u.deg) HEIGHT_ORM = u.Quantity(2199.835, u.m) location = EarthLocation.from_geodetic(lon=LON_ORM, lat=LAT_ORM, height=HEIGHT_ORM)